% Input

% *alpha* is the highest order of fractional derivative of given equation

% *beta* is the order of fractional derivatives other than alpha. *beta* must be a vector with decending ordered values

% *Uk* is the vector of coefficients

% *f unc* is defining the right hand side of given problem

% *t*0 and *R* denotes the left and right endpoints

% *y*0 is the initial conditions

% *m* denotes the number of steps

% *delta* is a real number greater than zero. We usually take *delta* = 1 or *delta* =fractional part of *alpha*

% Output

% *A* is an (*m* + 1) x (*m* + 1) matrix % *b* isan (*m* + 1) x 1matrix

% using *f ractionalTaylor*.*m*, where command *f ractionalTaylor*.*<sup>m</sup>* is defined by the Equation (18), gives us the linear system *AC* = *B* which is (*m* + 1)

% algebraic equations with unknown coefficients *C<sup>T</sup>*

% Next step is to use matlab function *linsolve*(*<sup>A</sup>*, *b*) to solve obtained algebraic equation for unknown coefficient vector *C<sup>T</sup>* with dimension (*m* + <sup>1</sup>). *C*= *linsolve*(*<sup>A</sup>*,*b*)

%Output

%*C* is an (*m* + 1) x 1 matrix which is the solution of linear system *AC* = *B*

% Next step is substituting obtained coefficients to *approxSoln*() as input, where the command *approxSoln*() defined by Equation (15), we ge<sup>t</sup> the approximate solution of given problem [*s*, *y*] = *approxSoln*(*C*) %

 Input

> % *C* is the vector of coefficients obtained in previous step.

% Output

> % *s* is the nodes on [*t*0, *R*] in which the approximate solution calculated % *y* is the numerical solution evaluated in the points of *s*.

#### **5. Illustrative Examples**

To illustrate the applicability and effectiveness of the presented method, we give six examples in this section. In each example, we apply the fractional Taylor operational matrix method which is presented in previous section and the approximate results compared with analytical solutions. Obtained results indicate that the proposed technique is very effective for multi-term fractional differential equations. In order to solve the numerical computations, MATLAB version R2015a has been used.

For choosing *δ*, we usually take either *δ* = 1 or *δ* = *α* − *<sup>α</sup>*, the fractional part of *α*.

## *5.1. Example 1*

Consider the following form of multi-order fractional differential equation [59]

$$\begin{aligned} D^a y(t) &= u\_0 D^{\bar{\mu}\_0} y(t) + u\_1 D^{\bar{\mu}\_1} y(t) + u\_2 D^{\bar{\mu}\_2} y(t) + u\_3 D^{\bar{\mu}\_3} y(t) + f(t), \ 0 \le t \le R, \\ y(0) &= V\_0, \ y'(0) = V\_1 \end{aligned} \tag{19}$$

We let *α* = 2, *V*0 = *V*1 = 0, *R* = 1, the coefficients *u*0 = *u*2 = −1, *u*1 = 2, *u*3 = 0 and *β*0 = 0, *β*1 = 1, *β*2 = 12 and the function *f*(*t*) is

$$f(t) = t^7 + \frac{2048}{429\sqrt{\pi}}t^{6.5} - 14t^6 + 42t^5 - t^2 - \frac{8}{3\sqrt{\pi}}t^{1.5} + 4t - 2.5$$

where the exact solution is *y*(*t*) = *t*7 − *t*2.

We apply the given procedure which is implemented in previous section for solving the Equation (19) step by step.

Firstly, change variable *t* ∈ [0, *R*] to *q* ∈ [0, 1] by using *q* = *t*/*R*.

Now, we use the Equation (8) and ge<sup>t</sup>

$$\frac{1}{R^a}D^a y\_1(q) = \frac{u\_0}{R^{\beta\_0}}D^{\beta\_0} y\_1(q) + \frac{u\_1}{R^{\beta\_1}}D^{\beta\_1} y\_1(q) + \frac{u\_2}{R^{\beta\_2}}D^{\beta\_2} y\_1(q) + \frac{u\_3}{R^{\beta\_3}}D^{\beta\_3} y\_1(q) + f\_1(q) \tag{20}$$

where 0 ≤ *q* ≤ 1. Next, using Equation (7) we ge<sup>t</sup>

$$\begin{split} \frac{1}{R^{\mathfrak{a}}} (y\_1(q) - y\_1(0) - q y\_1 \iota(0)) &= \frac{u\_0}{R^{\mathfrak{a}\_0}} I^{\mathfrak{a} - \beta\_0} (y\_1(q) - y\_1(0) - q y\_1 \iota(0)) \\ &+ \frac{u\_1}{R^{\mathfrak{a}\_1}} I^{\mathfrak{a} - \beta\_1} (y\_1(q) - y\_1(0) - q y\_1 \iota(0)) \\ &+ \frac{u\_2}{R^{\mathfrak{a}\_2}} I^{\mathfrak{a} - \beta\_2} (y\_1(q) - y\_1(0) - q y\_1 \iota(0)) \\ &+ \frac{u\_3}{R^{\mathfrak{a}\_3}} I^{\mathfrak{a} - \beta\_3} (y\_1(q) - y\_1(0) - q y\_1 \iota(0)) \\ &+ I^{\mathfrak{a}} f\_1(q). \end{split} \tag{21}$$

Now, using Equation (21) and substituting initial conditions *y*(0) = *V*0, *y*(0) = *V*1 into equation

$$\begin{split} \frac{1}{R^{\alpha}} (\mathbf{C}^{T} T\_{m\delta}(q) - V\_{0} - RqV\_{1}) &= \frac{u\_{0}}{R^{\beta\_{0}}} I^{\alpha - \beta\_{0}} (\mathbf{C}^{T} T\_{m\delta}(q) - V\_{0} - RqV\_{1}) \\ &+ \frac{u\_{1}}{R^{\beta\_{1}}} I^{\alpha - \beta\_{1}} (\mathbf{C}^{T} T\_{m\delta}(q) - V\_{0} - RqV\_{1}) \\ &+ \frac{u\_{2}}{R^{\beta\_{2}}} I^{\alpha - \beta\_{2}} (\mathbf{C}^{T} T\_{m\delta}(q) - V\_{0} - RqV\_{1}) \\ &+ \frac{u\_{3}}{R^{\beta\_{3}}} I^{\alpha - \beta\_{3}} (\mathbf{C}^{T} T\_{m\delta}(q) - V\_{0} - RqV\_{1}) \\ &+ I^{\alpha} f\_{1}(q). \end{split} \tag{22}$$

From Equation (12), we have

1 *Rα* (*CTTmδ*(*q*) − *V*0 − *RqV*1) = *u*0 *Rβ*0 *<sup>q</sup><sup>α</sup>*−*β*0*C<sup>T</sup>*(*<sup>G</sup><sup>α</sup>*−*β*<sup>0</sup> ∗ *Tmδ*(*q*)) − *<sup>u</sup>*0*q<sup>α</sup>*−*β*<sup>0</sup> *<sup>R</sup>β*0<sup>Γ</sup>(*α* − *β*0 + 1)*<sup>V</sup>*<sup>0</sup> − *<sup>u</sup>*0*q<sup>α</sup>*−*β*0+<sup>1</sup> *<sup>R</sup>β*0<sup>Γ</sup>(*α* − *β*0 + 2)*<sup>V</sup>*<sup>1</sup> + *u*1 *Rβ*<sup>1</sup> *<sup>q</sup><sup>α</sup>*−*β*1*C<sup>T</sup>*(*<sup>G</sup><sup>α</sup>*−*β*<sup>1</sup> ∗ *Tmδ*(*q*)) − *<sup>u</sup>*1*q<sup>α</sup>*−*β*<sup>1</sup> *<sup>R</sup>β*1<sup>Γ</sup>(*α* − *β*1 + 1)*<sup>V</sup>*<sup>0</sup> − *<sup>u</sup>*1*q<sup>α</sup>*−*β*1+<sup>1</sup> *<sup>R</sup>β*1Γ(*α* − *β*1 + 2)*<sup>V</sup>*<sup>1</sup> + *u*2 *Rβ*2 *<sup>q</sup><sup>α</sup>*−*β*2*C<sup>T</sup>*(*<sup>G</sup><sup>α</sup>*−*β*<sup>2</sup> ∗ *Tmδ*(*q*)) − *<sup>u</sup>*2*q<sup>α</sup>*−*β*<sup>2</sup> *<sup>R</sup>β*2<sup>Γ</sup>(*α* − *β*2 + 1)*<sup>V</sup>*<sup>0</sup> − *<sup>u</sup>*2*q<sup>α</sup>*−*β*2+<sup>1</sup> *<sup>R</sup>β*2<sup>Γ</sup>(*α* − *β*2 + 2)*<sup>V</sup>*<sup>1</sup> + *u*3 *Rβ*3 *<sup>q</sup><sup>α</sup>*−*β*3*C<sup>T</sup>*(*<sup>G</sup><sup>α</sup>*−*β*<sup>3</sup> ∗ *Tmδ*(*q*)) − *<sup>u</sup>*3*q<sup>α</sup>*−*β*<sup>3</sup> *<sup>R</sup>β*3<sup>Γ</sup>(*α* − *β*3 + 1)*<sup>V</sup>*<sup>0</sup> − *<sup>u</sup>*3*sα*−*β*3+<sup>1</sup> *<sup>R</sup>β*3<sup>Γ</sup>(*α* − *β*3 + 2)*<sup>V</sup>*<sup>1</sup> + *Iα f*1(*q*). (23)

Now, taking *R* = 1 in Equation (23) and putting the given values for *V*0, *V*1, *ui*, *βi* where *i* = 0, 1, 2, 3 into this equation, we ge<sup>t</sup>

$$\mathbf{C}^{T}T\_{m\delta} = 2q^{1}\mathbf{C}^{T}(\mathbf{G}\_{1} \ast T\_{m\delta}(q)) - q^{3/2}\mathbf{C}^{T}(\mathbf{G}\_{3/2} \ast T\_{m\delta}(q)) - q^{2}\mathbf{C}^{T}(\mathbf{G}\_{1} \ast T\_{m\delta}(q)) + l^{2}f\_{1}(q) \tag{24}$$

Finally, taking the collocation points *qj* = *j*/*m* (*j* = 0, 1, ..., *m*) generates a linear algebraic system of dimension *m* + 1 with unknown vector *CT*. In order to solve this system by using presented method and comparing the results, we choose *δ* = 1 and different values of *m*.

To show the efficiency, we compared the numerical results with the method given in [59].

Table 1, compares the obtained results for absolute error with *m* = 4, 6, 7. We observe from Table 1 that, the absolute errors for presented method are smaller and the numerical solution is more accurate for the same size of *m*.

**Table 1.** The comparison absolute errors of the present scheme and method given in [59] with *m* = 4, 6, 7.


In Figures 1–3, we present the graphical representation of comparison between exact solution and the numerical solutions obtained by proposed method and the method of [59] for the problem (19) with *m* = 4, 6, 7 respectively. From these results, we can conclude that *m* = 4 and *m* = 6 give larger absolute error, while *m* = 7 gives smaller absolute error (10−16) and more precise numerical solution. These comparisons also shows that the results obtained by proposed method is closer to the exact solution than the results of [59].

**Figure 1.** The comparison between exact solution and the numerical solutions obtained by proposed method and the method of [59] with *m*, *n* = 4.

**Figure 2.** The comparison between exact solution and the numerical solutions obtained by proposed method and the method of [59] with *m*, *n* = 6.

**Figure 3.** The comparison between exact solution and the numerical solutions obtained by proposed method and the method of [59] with *m*, *n* = 7.

In Figure 4, we show the graphical representation of absolute errors obtained by using proposed method and the method of [59] with *m*, *n* = 6.

**Figure 4.** The behaviour of absolute errors obtained by using proposed method and the method of [59] with *m*, *n* = 6.

From Figure 4, we can conclude that the absolute error obtained by our method is remaining smaller and stable while the absolute error of other method is increasing in the interval [0, 1].

In Figures 5 and 6, we give the graphical representation of absolute errors obtained by using proposed method with *m* = 4, 7 respectively.

**Figure 5.** The absolute error with *m* = 4.

**Figure 6.** The absolute error with *m* = 7.

A pseudo-code for MATLAB implementation of Example 1 is given in Algorithm 2 below :

**Algorithm 2:** Fractional Taylor Method

*alpha* = 2; *beta* = [1, 1/2, 0]; *Uk* = [2, −1, −<sup>1</sup>]; *f unc* =@(t) *t*7 + 2048/(429 ∗ *sqr<sup>t</sup>*(*p<sup>i</sup>*)) ∗ *t*6.5 − 14 ∗ *t*6 + 42 ∗ *t*5 − *t*2 − ... 8/(3 ∗ *sqr<sup>t</sup>*(*p<sup>i</sup>*)) ∗ *t*1.5 + 4 ∗ *t* − 2; *t*0 = 0 ; *R* = 1; *y*0 = [0; 0]; *m* = 4; *delta* = 1; [*<sup>A</sup>*, *b*] = *f ractionalTaylor*(*alpha*, *beta*, *Uk*, *f unc*, *t*0, *R*, *y*0, *m*, *delta*) *C* = *linsolve*(*<sup>A</sup>*, *b*) [*s*, *y*] = *approxSoln*(*C*)

*5.2. Example 2*

In this example, we consider the Equation (19) with *α* = 2, *V*0 = *V*1 = 0, the coefficients *u*0 = *u*2 = −1, *u*1 = 0, *u*3 = 2 and *β*0 = 0, *β*2 = 23 , *β*3 = 53 and the function is

$$f(t) = t^3 + 6t - \frac{12}{\Gamma(\frac{7}{3})}t^{\frac{4}{3}} + \frac{6}{\Gamma(\frac{10}{3})}t^{\frac{7}{3}}.$$

The exact solution of this equation is *y*(*t*) = *t*3 [59].

Applying the same procedure to given problem as presented in Example 1, we ge<sup>t</sup> the following equation

$$\mathbf{C}^{T}T\_{m\delta} = 2q^{1/3}\mathbf{C}^{T}(\mathbf{G}\_{1/3} \ast T\_{m\delta}(q)) - q^{4/3}\mathbf{C}^{T}(\mathbf{G}\_{4/3} \ast T\_{m\delta}(q)) - q^{2}\mathbf{C}^{T}(\mathbf{G}\_{2} \ast T\_{m\delta}(q)) + l^{2}f\_{1}(q) \tag{25}$$

As we stated in previous example, collocating this equation at the nodes *qj* = *j*/*m* (*j* = 0, 1, ..., *m*) generates a system of algebraic equations. In this example, to solve this sysem for *<sup>C</sup>T*, we choose *δ* = 1, 1.5 and different values of *m*.

Table 2 shows the results for obtained absolute errors by using presented method with *m* = 2, 3. From these results, we can see that, there is satisfactory agreemen<sup>t</sup> between the exact solution and numerical solutions. The absolute error is achieved about 10−15. We also note that, the proposed method gives better results for *m* = 2 by taking *δ* = 1.5.

**Table 2.** The absolute errors with *m* = 2, 3.


In Figure 7a, we show the graphical representation of obtained numerical solution and the exact solution of the given problem. Figure 7b presents the obtained absolute error by using proposed method with *m* = 3.

**Figure 7.** (**a**) The numerical and the exact solutions with *m* = 3. (**b**) The absolute error with *m* = 3.

#### *5.3. Example 3*

Consider the multi-term fractional order initial value problem [54]

$$D^{(2,2)}y(t) + 1.3D^{(1,5)}y(t) + 2.6y(t) = \sin(2t),\tag{26}$$

with initial conditions

$$y(0) = y(0) = y''(0) = 0,$$

where the equation have the series solution given by [52]

$$y\_s(t) = \frac{28561}{3600000}t^6 + \frac{2}{\Gamma(4.2)}t^{3.2} - \frac{13}{5\Gamma(4.9)}t^{3.9} + \frac{169}{50\Gamma(5.6)}t^{4.6}$$

$$-\frac{8}{\Gamma(6.2)}t^{5.2} - \frac{2197}{500\Gamma(6.3)}t^{5.3} - \frac{26}{5\Gamma(6.4)}t^{5.4} + \frac{52}{5\Gamma(6.9)}t^{5.9}.\tag{27}$$

In order to solve this problem, we choose *δ* = 1, and *m* = 10.

We give the comparison of series solution and the numerical solution obtained by presented method in Table 3. Table 4 compares the obtained absolute errors by using presented method with the results of [54]. From this compared results, it can be seen that the approximate solution is very close to series solution for a small number of *m* for the given method.

From the compared results of Table 4, we can conclude that the proposed method has better approach to series solution with a smaller *m*.


**Table 3.** Comparison of numerical solution with series solution for Example 3.

**Table 4.** Comparison of absolute errors for Example 3.


The graphical representation of comparison between series solution and numerical solutions obtained by presented method and the method of [54] in the interval [0, 1] is illustrated in Figure 8.

**Figure 8.** The comparison between series solution and numerical solutions obtained by proposed method and the method of [54] with *m* = 10.

In Figure 9, we show present graphical representation of absolute errors obtained by using proposed method and the method of [54] with *m* = 10.

**Figure 9.** The behaviour of absolute errors obtained by using proposed method and the method of [54].

In Figure 10, we show the graphical representation for series solution and the numerical results of presented method for the interval [0, <sup>10</sup>]. The results plotted in Figure 10 are in a very good and satisfactory agreemen<sup>t</sup> with the series solution given in [52] and the results of [60].

**Figure 10.** The behaviour of series solution and the numerical solution obtained by proposed method for the interval [0, <sup>10</sup>].

#### *5.4. Example 4*

Motivated by [50], we consider the following form of fractional differential equation,

$$D^a y(t) + y(t) = \begin{cases} \frac{2}{\Gamma(3-a)} t^{2-a} + t^2 - t, & a > 1\\ \frac{2}{\Gamma(3-a)} t^{2-a} - \frac{1}{\Gamma(2-a)} t^{1-a} + t^2 - t, & a \le 1 \end{cases} \tag{28}$$

with initial conditions

$$y(0) = 0, \ y(0) = -1$$

whose exact solution is *y*(*t*) = *t*2 − *t*.

In order to apply the presented method to Equation (28) and compare the results with methods of [54,61,62], we solve this problem with *α* = 0.3, 0.5, 0.7, 1.25, 1.5, 1.85, and different values for *δ* and *m*. The obtained results are presented as below.

In Table 5, we list the results of obtained absolute errors for *α* = 0.3, 0.5, 0.7 by use of presented method. Also, the results for *α* = 1.25, 1.5, 1.85 are given in Table 6.

**Table 5.** The absolute errors with *m* = 3 and *α* < 1 for Example 4.


**Table 6.** The absolute errors with *m* = 3 and *α* > 1 for Example 4.


In Figure 11a,b, we present the graphical representation of obtained results for numerical and exact solution of the given problem and absolute error for *α* = 1.5 in the interval [0, 1].

In Figure 12, we plot the graphical representation for behavior of the obtained numerical solution by use of the presented method and the exact solution of the given problem for *α* = 1.5 in the interval [0, <sup>15</sup>].

**Figure 11.** (**a**) The numerical and exact solutions for *α* = 1.5. (**b**) The absolute error for *α* = 1.5.

**Figure 12.** The behaviour of the obtained numerical and exact solutions with *α* = 1.5 for the interval *t* ∈ [0, <sup>15</sup>].

Table 7 lists the obtained absolute errors for the given problem (28) at *t* = 1, 5, 10, 50 and *α* = 1.5 by use of presented method and some other methods in literature [54,61,62]. From this compared results, we can say that the numerical solution obtained by use of proposed method is in better agreemen<sup>t</sup> with the exact solution and obtained absolute error is smaller.

**Table 7.** Comparison of absolute errors between proposed method and some other numerical methods in literature at *t* = 1, 5, 10, 50 for *α* = 1.5.


In Figure 13, the behaviour of absolute error for *α* = 1.5 with *m* = 4 and *δ* = 1/2, 1 at *t* ∈ [0, 50] is presented. From this graph, it can be seen that we ge<sup>t</sup> better results by taking *δ* = 1/2 for this example and the numerical solution is very close to exact solution for a small number of *m*.

**Figure 13.** The behaviour of the absolute errors for proposed method where *α* = 1.5, *t* ∈ [0, 50] with *m* = 4 and *δ* = 1/2, 1.

### *5.5. Example 5*

In this example, we consider the following form of linear multi-term fractional differential equation with variable coefficients [65]

$$aD^2y(t) + b(t)D^{\beta\_1}y(t) + c(t)Dy(t) + e(t)D^{\beta\_2}y(t) + k(t)y(t) = f(t),\tag{29}$$

with,

$$y(0) = 2, \ y(0) = 0$$

where 0 < *β*2 < 1, 1 < *β*1 < 2 and

$$f(t) = -a - \frac{b(t)}{\Gamma(3-\beta\_1)}t^{2-\beta\_1} - c(t)t - \frac{c(t)}{\Gamma(3-\beta\_2)}t^{2-\beta\_2} + k(t)\left(2 - \frac{t^2}{2}\right)t$$

whose the exact solution is given by *y*(*t*) = 2 − *t*22

We give the numerical solution for the given problem by proposed method for *a* = 1, *b*(*t*) = √*<sup>t</sup>*, *c*(*t*) = *t* 13 ,*<sup>e</sup>*(*t*) = *t* 14 , *k*(*t*) = *t* 15 , *β*2 = 0.333, *β*1 = 1.234 with *δ* = 1.

 .

In Table 8, we give the results for maximum errors obtained by use of proposed method and comparison with the results of [65,66]. From this compared results, we can see that the numerical solution obtained by use of proposed method is closer to the exact solution.

**Table 8.** Maximum errors of Example 5 for *R* = 1 with *m* = 3, 4, 5, 6, 10, 20, 40.


Figure 14 presents the graphical representation for behaviour of numerical and exact solutions with *m* = 6. From this representation, we can see that the numerical solution is in a very good agreemen<sup>t</sup> with exact solution.

**Figure 14.** The behaviour of the numerical and exact solutions with *m* = 6.

#### *5.6. Example 6*

For the last example, let us consider the below fractional differential equation [63]

$$y(t) + D^{1/2}y(t) - 2y(t) = 0, \; t \in (0, R], \tag{30}$$

$$y(0) = 1$$

which arises, for example, in the study of generalized Basset force occuring when a spherical object sinks in a (relatively dense) incompressible viscous fluid; see [43,67]. By use of Laplace transformation of Caputo derivatives, we ge<sup>t</sup> the analytical solution as following

$$y(t) = \frac{2}{3\sqrt{t}} E\_{1/2, 1/2}(\sqrt{t}) - \frac{1}{6\sqrt{t}} E\_{1/2, 1/2}(-2\sqrt{t}) - \frac{1}{2\sqrt{\pi t}} e$$

where the Mittag–Leffler function *<sup>E</sup>λ*,*<sup>μ</sup>*(*t*) with parameters *λ*, *μ* > 0 is given as

$$E\_{\lambda,\mu}(t) = \sum\_{k=0}^{\infty} \frac{t^k}{\Gamma(\lambda k + \mu)}.$$

This Mittag–Leffler function and its variations are very significant in fractional calculus and fractional differential equations [68].

In order to solve given problem by use of proposed method and compare the results, we take *t* ∈ (0, 5] and use different values of *δ* and *m*.

Table 9 lists the exact and obtained numerical solutions by use of presented method and method of [63] for the given problem for *m* = 5, 10, 15, 20. Comparison of this results shows that, even for small values of *m*, the numerical solution obtained by use of presented method is in a better agreemen<sup>t</sup> with exact solution.


**Table 9.** The resulting values for Example 6, with *R* = 5 in some values of *t.*

In Figures 15a, 16a and 17a, we present the graphical representation of comparison between exact solution and the numerical solutions obtained by using proposed method and the method of [63] with taking *m* = 5, 10, 20 respectively. Also in Figures 15b, 16b and 17b we show the behaviour of absolute errors obtained by proposed method and the method of [63] in the interval [0, 1] with *m* = 5, 10, 20.

**Figure 15.** (**a**) The comparison of analytical solution and numerical solutions obtained by the proposed method and the method of [63] with *m* = 5. (**b**) The behaviour of the absolute errors between the exact solution and numerical solutions obtained by our method and the method given in [63] with *m* = 5.

**Figure 16.** (**a**) The comparison of analytical solution and numerical solutions obtained by the proposed method and the method of [63] with *m* = 10. (**b**) The behaviour of the absolute errors between the exact solution and numerical solutions obtained by our method and the method given in [63] with *m* = 10.

**Figure 17.** (**a**) The comparison of analytical solution and numerical solutions obtained by the proposed method and the method of [63] with *m* = 20. (**b**) The behaviour of the absolute errors between the exact solution and numerical solutions obtained by our method and the method given in [63] with *m* = 20.

From these graphical results represented in Figures 15–17, we can conclude that the absolute error obtained by our method is remaining smaller when compared the absolute error of method given in Reference [63].

## **6. Conclusions**

In this work, an operational matrix based on the fractional Taylor vector is used to numerically solve the multi-term fractional differential equations by reducing them to a set of linear algebraic equations, which simplifies the problem. From comparison of the obtained results with exact solutions and also with results of other methods in the literature, we conclude that the proposed method provides the solution with high accuracy. The findings also show that, even for the small number of steps, we can ge<sup>t</sup> satisfactory results by using presented method. All computational results are obtained by using MATLAB.

**Author Contributions:** Formal analysis, ˙ I.A.; Supervision, N.I.M. All authors contributed equally to this article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.
