**1. Introduction**

Fractional calculus is an emerging field of mathematics, which is a generalisation of differentiation and integration to non-integer orders. The history of fractional calculus is almost as long as the history of classical calculus, beginning with some speculations of Leibniz (1695, 1697) and Euler (1730). However, fractional calculus and fractional differential equations (FDEs) are increasingly becoming popular in recent years. The progressively developing history of this old and ye<sup>t</sup> novel topic can be found in [1–5]. In fact, fractional calculus provides the mathematical modeling of some important phenomena like social and natural in a more powerful way than the classical calculus. During the last few decades, many applications were reported in many branches of science and engineering such as chaotic systems [6,7], fluid mechanics [8], viscoelasticity [9], optimal control problems [10,11], chemical kinetics [12,13], electrochemistry [14], biology [15], physics [16], bioengineering [17], finance [18], social sciences [19], economics [20,21], optics [22], chemical reactions [23], rheology [24], and so on. Due to the importance of FDEs, the solutions of them are attracting widespread interest. On the other hand, analytical solutions are not always possible for solving them. Therefore, numerical techniques becomes more important for solving such equations.

There are various numerical methods have been developed for solving FDEs in literature such as predictor-corrector method [25], Laplace transforms [26], Taylor collocation method [27], variational iteration method and homotopy perturbation method [8] (Chapter 6), Adomian decomposition method [28], Tau method [29], inverse Laplace transform [30], Haar wavelet collocation method [31], generalized block pulse operational matrix [32], shifted Legendre-tau method [33], fractional multi-step differential transformed method [34], q-homotopy analysis transform method [35], conformable Laplace transform [36], fractional B-splines collocation method [37], finite difference method [38], homotopy analysis method [39] and so on.

Multi-term fractional differential equations are one of the most important type of FDEs, which is a system of mixed fractional and ordinary differential equations and involving more than one fractional differential operators. Nowadays, they are widely appearing for modelling of many important processes, especially for multirate systems. Their numerical solution is then a strong subject that deserves high attention. In this paper, motivated by the results reported in [40,41] for solving a smaller class of problems where the highest order of derivative is an integer and involving at most one noninteger order derivative, we go further and establish a method for numerical solutions for higher order and arbitrary multi-term fractional differential equations which have a general form

$$D^{\mathfrak{a}}y(t) = f\left(t, y(t), D^{\mathfrak{f}\_0}y(t), D^{\mathfrak{f}\_1}y(t), \dots, D^{\mathfrak{f}\_k}y(t)\right), \ t \in [0, R] \tag{1}$$

where *Dα* representing the Caputo fractional derivative of order *α* > 0 and we assume that 0 < *β*0 < *β*1 < ... < *βk* < *α*, *y*(*p*) = *Yp*, *p* = 0, 1, ...*n* where *n* − 1 < *α* < *n*.

Multi-term fractional order differential equations also have useful properties and they can describe complex multi-rate physical processes in a various way and can be applied in many fields, see e.g., [2,4,26,42]. Basset [43] and Bagley–Torvik [44] equations can be given as important examples for smaller class of multi-term fractional differential equations. Existence, uniqueness and stability of solution for multi-term fractional differential equations are discussed in [45–49]. Because of difficulty of finding the exact solutions for such equations, many new numerical techniques have been developed to investigate the numerical solutions such as Adams method [50], Haar wavelet method [51], differential transform method [52], Adams–Bashforth–Moulton method [53], collocation method based on shifted Chebyshev polynomials of the first kind [54], Boubaker polynomials method [55], matrix Mittag–Leffler functions [56], differential transform method [57] and so on.

Our main purpose is to present an effective, reliable method to approximate initial value problem for the Equation (1). In order to reach this aim, we rewrite and focus the general type of Caputo multi-term fractional differential equation given in Equation (1) in the following linear form

$$D^\pi y(t) = \sum\_{i=0}^k u\_i D^{\beta\_i} y(t) + u\_{k+1} y(t) + f(t), \ 0 \le t \le R,\tag{2}$$

subject to the

$$\begin{aligned} y^{(p)}(0) &= \Upsilon\_{p\prime} \ p = 0, 1, \ldots, n-1 \text{ where } n-1 < \alpha < n\\ u\_i \ (i = 0, 1, \ldots, k) \text{ are known coefficients and} \\ 0 &< \beta\_0 < \beta\_1 < \ldots < \beta\_k < n \end{aligned} \tag{3}$$

Here, we also state that the highest order *α* need not to be an integer. This equation is important in applications due to the fact it can treat the problems with fractional force, therefore it is suitable for being treated within fractional operators of Caputo type.

In this work, a numerical approach based on fractional Taylor vector is proposed to solve the initial value problem of general type of multi-term fractional differential equations which is given in Equations (2) and (3). The core idea of this method is to employ the operational matrix of fractional integration based on fractional Taylor vector to given problem and reduce it to a set of algebraic equations which can be efficiently solved.

The structure of the manuscript is organized as follows. In Section 2, we briefly introduce some preliminary ideas of fractional calculus and necessary definitions. In Section 3, an operational matrix of fractional integration based on fractional taylor vector is derived. In Section 4, we present the numerical algorithm to solve the given equation and a pseudo-code for matlab is also provided in Algorithm 1. In Section 5, the presented method is applied to six examples to demonstrate the efficiency. A final conclusion is presented in the last section.

#### **2. Preliminary Knowledge**

In this section, we recall some fundamental definitions and preliminary facts of fractional calculus.

#### *2.1. The Fractional Integral and Derivative*

**Definition 1.** *The Riemann–Liouville fractional integral to order α of an integrable function y*(*t*) *is defined to be*

$$M^a y(t) = \begin{cases} \frac{1}{\Gamma(a)} \int\_0^t (t-s)^{a-1} y(q) \, \text{ds}, & a > 0\\ y(t), & a = 0 \end{cases} \tag{4}$$

*When applied to a power function, it yields the following result:*

$$I^a(t)^c = \frac{\Gamma(c+1)}{\Gamma(c+a+1)} (t)^{c+a}, \ a \ge 0, \ c > -1\tag{5}$$

*The operator has a semigroup property, namely*

$$I^{\alpha}I^{\beta}y(t) = I^{\beta}I^{\alpha}y(t), \; \alpha, \beta > 0$$

*and it is linear, namely*

$$I^a(A\_1 y\_1(t) + A\_2 y\_2(t)) = A\_1 I^a y\_1(t) + A\_2 I^a y\_2(t)$$

*for any two functions y*1*,y*2 *and constants A*1*,A*2*.*

**Definition 2.** *The fractional derivative of y*(*t*) *of the order α in the Caputo sense is given as*

$$D^a y(t) = l^{j-a} \left(\frac{\mathbf{d}^j}{\mathbf{d}t^j} y(t)\right), \ j - 1 < a \le j, \ j \in \mathbb{N} \tag{6}$$

#### *2.2. Some Properties*

1. The Riemann-Liouville fractional integral and Caputo fractional derivative do not usually commute with each other. The following Newton–Leibniz identity gives an important relation between them:

$$d^a(D^a y(t)) = y(t) - \sum\_{i=0}^{j-1} y^{(i)}(0) \frac{t^i}{i!} \tag{7}$$

2. The Caputo fractional derivative also has the following substitution identity. If we write *y*1(*q*) = *y*(*q<sup>R</sup>*) and *q* = *t*/*R*, then

$$D^a y(t) = \frac{1}{R^a} D^a y\_1(q) \tag{8}$$

where *j* − 1 < *α* ≤ *j*, *j* ∈ N

#### **3. Operational Matrix of Fractional Integration for Fractional Taylor Vector**

#### *3.1. Fractional Taylor Basis Vector*

We shall make use of the fractional Taylor vector,

$$T\_{m\delta}(t) = \begin{bmatrix} 1, t^{\delta}, t^{2\delta}, \dots, t^{m\delta} \end{bmatrix} \tag{9}$$

for *m* ∈ N and *δ* > 0 in the work of this paper.

#### *3.2. Approximation of Function*

Suppose that *Tmδ*(*t*) ⊂ *H*, where *H* is the space of all square integrable functions on the interval [0, 1]. For any *y* ∈ *H*, since *S* = *span* /1, *tδ*, *<sup>t</sup>*2*δ*, ..., *t<sup>m</sup><sup>δ</sup>*0 is a finite dimensional vector space in *H*, then, *y* has a unique best approximation *y*∗ ∈ *S*, so that

$$\forall \hat{y} \in \mathcal{S}\_{\prime} \; ||y - y\_{\ast}|| \le ||y - \hat{y}||$$

Therefore, the function *y* is approximated by fractional Taylor vector as following

$$y \simeq y\_\* = \sum\_{i=0}^{m} c\_i t^{i\delta} = \mathbb{C}^T T\_{m\delta}(t) \tag{10}$$

where *Tmδ*(*t*) denote the fractional Taylor vector and

$$\mathbf{C}^{T} = [\mathbf{c}\_{0}, \mathbf{c}\_{1}, \mathbf{c}\_{2}, \dots, \mathbf{c}\_{m}] \tag{11}$$

are the unique coefficients.

#### *3.3. Fractional Taylor Operational Matrix of Integration*

By using the property of Riemann-Liouville fractional integral given in Equations (5) and (9), we ge<sup>t</sup>

$$\begin{split} I^{\mathfrak{a}}(T\_{m\delta}(t)) &= \left[ \begin{array}{c} \frac{1}{\Gamma(a+1)} t^{\mathfrak{a}}, \frac{\Gamma(\delta+1)}{\Gamma(\delta+a+1)} t^{\delta+a}, \frac{\Gamma(2\delta+1)}{\Gamma(2\delta+a+1)} t^{2\delta+a}, \dots, \frac{\Gamma(m\delta+1)}{\Gamma(m\delta+a+1)} t^{m\delta+a} \right] \\ &= t^{\mathfrak{a}} M\_{\mathfrak{a}} T\_{m\delta}(t) \end{split} \tag{12}$$

where

$$M\_{\mathfrak{a}} = \operatorname{diag} \left[ \frac{1}{\Gamma(\mathfrak{a}+1)}, \frac{\Gamma(\delta+1)}{\Gamma(\delta+\mathfrak{a}+1)}, \frac{\Gamma(2\delta+1)}{\Gamma(2\delta+\mathfrak{a}+1)}, \dots, \frac{\Gamma(m\delta+1)}{\Gamma(m\delta+\mathfrak{a}+1)} \right]$$

denotes the operational matrix of integration.

> If we define *Gα* as

$$G\_{\mathfrak{a}} = \left[ \frac{1}{\Gamma(\mathfrak{a}+1)}, \frac{\Gamma(\delta+1)}{\Gamma(\delta+\mathfrak{a}+1)}, \frac{\Gamma(2\delta+1)}{\Gamma(2\delta+\mathfrak{a}+1)}, \dots, \frac{\Gamma(m\delta+1)}{\Gamma(m\delta+\mathfrak{a}+1)} \right]$$

then, we can rewrite the Equation (10) as

$$I^{\mathfrak{a}}(T\_{m\delta}(t)) = t^{\mathfrak{a}} G\_{\mathfrak{a}} \ast T\_{m\delta}(t) \tag{13}$$

where ∗ denotes the operation of multiplying matrices term by term.

#### **4. The Numerical Algorithm**

In this section, to solve the given multi-term fractional differential equation in Equations (2) and (3), we employ the fractional Taylor method. The algorithm of method is given below.

Firstly, by using the transformation *q* = *t*/*R*, we replace the variable *t* ∈ [0, *R*] with *q* ∈ [0, 1]. Now, by using Equation (8) in Equation (2), we ge<sup>t</sup>

$$\frac{1}{R^a}D^a y\_1(q) = \sum\_{i=0}^k \frac{1}{R^{\beta\_i}} \mu\_i D^{\beta\_i} y\_1(q) + \mu\_{k+1} y\_1(q) + f\_1(q), 0 \le s \le 1 \tag{14}$$

where *f*1(*q*) = *f*(*q<sup>R</sup>*) and *y*1(*q*) = *y*(*q<sup>R</sup>*). Similar to Equation (10) we approximate the *y*1(*q*) as

$$y\_1(q) = \sum\_{i=0}^{m} c\_i q^{i\delta} = \mathbb{C}^T T\_{m\delta}(q) \tag{15}$$

such that *Tmδ*(*q*)=[1, *qδ*, *<sup>q</sup>*<sup>2</sup>*δ*, ..., *q<sup>m</sup><sup>δ</sup>*]*<sup>T</sup>* is the fractional Taylor vector and the unique coefficients *C<sup>T</sup>* is given in Equation (11).

Next, applying the Riemann–Liouville fractional integral on both side of (14), we ge<sup>t</sup>

$$\frac{1}{R^a} \left[ y\_1(q) - \sum\_{j=0}^{n-1} y\_1^{(j)}(0^+) \frac{t^j}{j!} \right] = \sum\_{i=0}^k \frac{1}{R^{\beta\_i}} u\_i I^{a-\beta\_i} \left[ y\_1(q) - \sum\_{j=0}^{n\_i-1} y\_1^{(j)}(0^+) \frac{t^j}{j!} \right]$$

$$+ u\_{k+1} I^a y\_1(q) + I^a f\_1(q) \tag{16}$$

where *y*(*p*)(0) = *Vp*, *p* = 0, 1, ..., *n* − 1 where *ni* − 1 < *βi* < *ni*.

Hence, by substituting initial conditions (3), we ge<sup>t</sup>

$$\frac{1}{R^a} \left[ y\_1(q) \right] = \sum\_{i=0}^k \frac{1}{R^{\beta\_i}} u\_i I^{a-\beta\_i} \left[ y\_1(q) \right] + u\_{k+1} I^a y\_1(q) + h\_1(q) \tag{17}$$

such that *h*1(*q*) = *Iα f*1(*q*) + 1*Rα* ∑*<sup>n</sup>*−<sup>1</sup> *j*=0 *Vj tjj*! + ∑*ki*=<sup>0</sup> 1*Rβi uiI<sup>α</sup>*−*β<sup>i</sup>* ∑*ni*−<sup>1</sup> *j*=0 *Vj tjj*! .

Now, by using the Equation (12), we approximate the fractional order integrals in Equation (17) and we have

$$\frac{1}{R^{\alpha}} \left[ \mathbb{C}^{T} T\_{m\delta}(q) \right] = \sum\_{i=0}^{k} \frac{1}{R^{\beta\_i}} \mu\_i \mathbb{C}^{T} q^{\alpha - \beta\_i} \left( \mathbb{G}\_{a - \beta\_i} \* T\_{m\delta}(q) \right)$$

$$+ \mu\_{k+1} q^a \mathbb{C}^{T} (\mathbb{G}\_{a} \* T\_{m\delta}(q)) + h\_1(q) \tag{18}$$

Finally, by taking the collocation points *qj* = *j*/*m* (*j* = 0, 1, ..., *m*) in Equation (18), we ge<sup>t</sup> *m* + 1 linear algebraic equations. This linear system can be solved for the unknown vector *CT*. Consequently, *y*1(*q*) can be approximated by Equation (15).

#### *MATLAB Implementation of Method*

The pseudocode given in Algorithm 1 below allows us to use proposed method in MATLAB for obtain a numerical solution of given problem [58].

**Algorithm 1:** Fractional Taylor Method

[*A*, *b*] = *f ractionalTaylor*(*alpha*, *beta*, *Uk*, *f unc*, *t*0, *R*, *y*0, *m*, *delta*)
