During the past decades, the field of fractional differential equations attracted the interest of researchers in several areas including physics, chemistry, engineering and even finance and social sciences (see [1, 2]) and there was significant interest in developing numerical schemes for their solution. These methods include Laplace transforms (see [3]), Fourier transforms (see [4]), eigenvector expansion (see [5]), adomian decomposition method (ADM) (see [6, 7]), variation iteration method (VIM)(see [8, 9]), fractional differential transform method (FDTM) (see [10, 11]), fractional difference method (FDM) (see [12]) and power series method (see [13]). But, few papers were reported to solve the fractional order differential equations by application of wavelets (see [14, 15]).

We intend to use the Chebyshev wavelet method to solve multi-order arbitrary differential equations with variable coefficients. First, we construct the Chebyshev wavelets and derive the operational matrix of fractional integration; then the underlying fractional differential equation is converted to an algebraic equation by introducing the wavelet operational matrix of the fractional integration. In this paper, by using the Chebyshev wavelets, we solve the following multi-order fractional differential equations (FDEs) with variable coefficients

$\begin{eqnarray}\label{eq:1)} &&{D_{*}^{v}}u(t) + \sum\limits_{i = 1}^{r - 1} {{\gamma_i}(t){D_ * }^{{\beta _i}}u(t)} + {\gamma _r}(t)u(t)\ = g(t)\quad, t\in[0, 1), \end{eqnarray}$

(1.1)

$\begin{eqnarray}\label{eq:2)} &&{u^{(i)}}( 0 )={d_i}, \quad i= 0, 1, \cdots, n-1, \end{eqnarray}$

(1.2)

where$\quad 0 < {\beta _1} < {\beta _2} < \cdots <{\beta_{r-1}} <v, \quad n-1 <v\leq n$ are constants. Moreover, ${D_{*}^{v}}$ denotes the Caputo fractional derivative of order ${v}$, the values of $ {d_{i}}{(i=0, 1 \cdots, n-1)}$describe the initial state of ${u(t)}$ and ${g(t)}$ is a known function. The existence and uniqueness of solutions of FDEs were studied by [16].

The paper is organized as follows. In Section 2, we introduce some necessary definitions and mathematical preliminaries of fractional calculus. In Section 3, after describing the basic formulation of wavelets and the Chebyshev wavelets, we derive the Chebyshev wavelet operational matrix of the fractional differential equation. In Section 4, the method is defined for approximate solution of the fractional problem (1.1) and (1.2). In Section 5, the error analysis technique based on the residual function is developed for the present method. A numerical example is given to demonstrate the validity of the method in solving fractional differential equation in Section 6. Section 7 comments on the result.

We give some necessary deflnitions and mathematical preliminaries of the fractional calculus theory which are used further in this paper.

Definition 2.1  The Riemann-Liouville fractional integral operator ${I^\alpha}$ of order $\alpha(\alpha>0)$ on usual Lebesgue space ${L^2{[a, b]}}$ is given by

$\begin{eqnarray}\label{eq:3)} {I^\alpha }f(t) = \frac{1}{{\Gamma (\alpha )}}\int_{\, a}^{\, t} {{{(t - \tau )}^{\alpha - 1}}f(\tau )d\tau, t > a}, \end{eqnarray}$

(2.1)

$\begin{eqnarray}\label{eq:4)} {I^0}f(t) = f(t), \end{eqnarray}$

(2.2)

and its fractional derivative of order $\alpha>0$ is normally used:

$\begin{equation}\label{eq:5)} {D^\alpha }f(t) = \frac{{{d^n}}}{{d{t^n}}}({I^{n -\alpha }}f(t)), \quad n - 1 < \alpha \leq n, \end{equation}$

(2.3)

where $n$ is an integer. For Riemann-Liouville definition, one has

$\begin{equation}\label{eq:6)} {I^\alpha }{t^\nu } = \frac{{\Gamma (\nu +1)}}{{\Gamma (\alpha + \nu + 1)}}{t^{\alpha + \nu }}. \end{equation}$

(2.4)

The Riemann-Liouville derivative has certain disadvantages when trying to model real-world phenomena with fractional differential equations. Therefore, we shall introduce now a modified fractional differential operator $D_{*}^{\alpha}$ proposed by Caputo.

Definition 2.2  The Caputo definition of fractional differential operator is given by

$\begin{equation}\label{eq:7)} {D_ * }^\alpha f(t) = \frac{1}{{\Gamma (n -\alpha )}}\int_{\, a}^{\, t} {{{(t - \tau )}^{n - \alpha - 1}}{f^{(n)}}(\tau )d\tau, }\quad n -1 <\alpha \leq n, \end{equation}$

(2.5)

where $n$ is an integer.

It has the following basic properties for $n -1 <\alpha \leq n$ and $f\in{L^2{[a, b]}}$,

$\begin{eqnarray}\label{eq:8)} &&{D_ * }^\alpha {I^\alpha }f(t) = f(t), \end{eqnarray}$

(2.6)

$\begin{eqnarray}\label{eq:9)} &&\begin{aligned} {I^\alpha }{D_ * }^\alpha f(t) = f(t) - \sum\limits_{k= 0}^{n - 1} {{f^{(k)}}({a^ + })\frac{{{t^k}}}{{k!}}}, t> a, \end{aligned}\end{eqnarray}$

(2.7)

$\begin{eqnarray}\label{eq:10)} &&D_ * ^\beta f(t) = {I^{v - \beta }}D_ * ^vf(t). \end{eqnarray}$

(2.8)

In this section, we use the Chebyshev polynomials to construct the Chebyshev wavelets and give some properties of the wavelets.

Wavelets are a family of functions constructed from dilation and translation of a single function $\psi(t)$ called the mother wavelet. When the dilation parameter $a$ and the translation parameter $b$ vary continuously we have the following family of continuous wavelets as [17]:

$\begin{equation}\label{eq:11)} {\psi _{ab}}(t) = {\left| a \right|^{ - \frac{1}{2}}}\psi (\frac{{t - b}}{a}), \;a, b \in R, a \neq 0. \end{equation}$

(3.1)

If we restrict the parameters $a$ and $b$ to discrete values as $a = a_0^{ - k}, b = n{b_0}a_0^{ - k}$, ${a_0} > 0$, $b_0 > 0$, where $n$ and $k$ are positive integers, the family of discrete wavelets are defined as

$\begin{equation}\label{eq:12)} {\psi _{kn}}(t) = {\left| {{a_0}}\right|^{\frac{k}{2}}}\psi (a_0^kt - n{b_0}), \end{equation}$

(3.2)

where $\psi_{kn}$ form wavelet bases for $L^2{(R)}$. In particular, when $a_0 = 2\;\text{and}\;b_0 = 1$, $\psi_{kn}$ form orthogonal bases.

The Chebyshev wavelets ${\psi _{nm}}(t) = \psi (k, n, m, t)$, which are defined on the interval $[0, 1)$, involve four arguments. That is, ${k}$ is assumed any positive integer, $n = 1, 2, \cdots, {2^{k- 1}}$, $m$ is the degree of the Chebyshev polynomials and $t$ is the normalized time. They are defined on the interval $[0, 1)$ as

$\begin{equation}\label{eq:13)} \psi _{nm}(t) =\left\{ \begin{aligned} {2^{\frac{k}{2}}}{{\tilde T}_m}({2^k}t- 2n + 1), &\frac{{n- 1}}{{{2^{k- 1}}}} \leq t <\frac{n}{{{2^{k - 1}}}}, \\ 0, \quad \quad \quad \quad \quad \quad \quad \quad &\text{otherwise, }\\ \end{aligned} \right. \end{equation}$

(3.3)

where

$\begin{equation}\label{eq:14)} \begin{aligned} &{\tilde T_0} =\frac{1}{{\sqrt \pi }}, \\ &{\tilde T_m}(t) = \sqrt {\frac{2}{\pi }} {T_m}(t), m = 1, 2, \cdots, M - 1, \end{aligned} \end{equation}$

(3.4)

here $T_m(t)$ are the Chebyshev polynomials of degree $m$ which respect to the weight function $\omega (t) = 1/\sqrt {1 - {t^2}}$ on interval [-1,1] and satisfy the following recursive formula

$\begin{equation}\label{eq:15)} {\begin{aligned} T_0(t) = 1, T_{1}(t)=t, T_{{m+1}}(t)=2t{T_m}(t)-{T_{m-1}}(t), m=1, 2, \cdots. \end{aligned}} \end{equation}$

(3.5)

A function $f(t)$ defined on the interval [0,1] may be expanded as

$\begin{equation}\label{eq:16)} f(t) = \sum\limits_{n = 1}^\infty {\sum\limits_{m = 0}^\infty {{c_{nm}}{\psi _{nm}}(t)} }, \end{equation}$

(3.6)

where

$\begin{equation}\label{eq:17)} {c_{nm}} = {(f(t), {\psi _{nm}}(t))_{{\omega _n}}} = \int_0^1 {{\omega _n}(t){\psi _{nm}}(t)f(t)dt}. \end{equation}$

(3.7)

If the infinite series in eq. (3.6) is truncated, then eq. (3.6) can be written as

$\begin{equation}\label{eq:18)} f(t) \approx \sum\limits_{n = 1}^{{2^{k - 1}}} {\sum\limits_{m = 0}^{M - 1} {{c_{nm}}{\psi _{nm}}(t)} } = {C^T}\Psi (t), \end{equation}$

(3.8)

where $T$ indicates transposition, $C$ and $\Psi(t)$ are $2^{k-1}M$ dimensional column vectors given by

$\begin{equation}\label{eq:19)} C = {[{c_{10}}, {c_{11}}, \cdots, {c_{1(M-1)}}, {c_{20}}, \cdots, {c_{2(M-1)}}, \cdots, {c_{{2^{k-1}}0}}, \cdots, {c_{{2^{k - 1}}(M - 1)}}]^T}, \end{equation}$

(3.9)

$\begin{equation}\label{eq:20)} \Psi (t) = {[{\psi _{10}}(t), \cdots, {\psi _{1(M-1)}}(t), {\psi _{20}}(t), \cdots, {\psi _{2(M-1)}}(t), \cdots, {\psi _{{2^{k-1}}0}}(t), \cdots, {\psi _{{2^{k - 1}}(M - 1)}}(t)]^T}. \end{equation}$

(3.10)

Taking the collocation points as following

$\begin{equation}\label{eq:21)} {t_i} = \frac{{2i - 1}}{{{2^k}M}}, i = 1, 2, \cdots, {2^{k - 1}}M. \end{equation}$

(3.11)

We define the Chebyshev wavelet matrix ${\Phi _{m' \times m'}}$ as

$\begin{equation}\label{eq:22)} {\Phi _{m' \times m'}} = [\Psi ({t_1}), \Psi ({t_2}), \cdots, \Psi ({t_{m'}})], \end{equation}$

(3.12)

where $m' = {2^{k - 1}}M$.

The integration of the vector $\Psi(t)$ defined in eq. (3.10) can be obtained as

$\begin{equation}\label{eq:23)} \int_0^1 {\Psi (t)dt} \approx {\rm P}\Psi (t), \end{equation}$

(3.13)

where $P$ is the ${2^{k - 1}M \times {2^{k - 1}}}M$ operational matrix for integration [17].

Our purpose is to derive the Chebyshev wavelet operational matrix of the fractional integration. For this purpose, we rewrite Riemann-Liouville fractional integration, as following

$\begin{equation}\label{eq:24)} {I^\alpha }f(t) = \frac{1}{{\Gamma (\alpha )}}\int_{\, 0}^{\, t} {{{(t - \tau )}^{\alpha - 1}}f(\tau )d\tau = \frac{1}{{\Gamma (\alpha )}}{t^{\alpha - 1}} * f(t), t > 0}. \end{equation}$

(3.14)

Now, if $f(t)$ is expanded in the Chebyshev wavelets, as showed in eq. (3.8), the Riemann-Liouville fractional integration becomes

$\begin{equation}\label{eq:25)} {I^\alpha }f(t) = \frac{1}{{\Gamma (\alpha )}}{t^{\alpha - 1}} * f(t) \approx {C^T}\frac{1}{{\Gamma (\alpha )}}\left\{ {{t^{\alpha - 1}} * \Psi (t)} \right\}. \end{equation}$

(3.15)

Thus if ${t^{\alpha - 1}} * f(t)$ can be integrated, then expanded in the Chebyshev wavelets, the Riemann-Liouville fractional integration is solved.

Also, we define a $m'$-set of Block Pulse function (BPF) as

$\begin{equation}\label{eq:26)} b _i(t) =\left\{ \begin{aligned} 1, \quad &\frac{i}{m'} \leq t <\frac{i+1}{m'}, \\ 0, \quad &\text{otherwise}, \end{aligned} \right. \end{equation}$

(3.16)

where $i=0, 1, 2, \cdots, m'-1$.

The functions $b_i(t)$ are disjoint and orthogonal. That is

$\begin{eqnarray}\label{eq:27)} &&b_i(t)b_l(t)=\left\{ \begin{aligned} b_i(t), \quad i=l, \\ 0, \quad \quad i\neq l, \\ \end{aligned} \right.\end{eqnarray}$

(3.17)

$\begin{eqnarray}\label{eq:28)} &&\int_0^1 {}b_i(t)b_l(t)=\left\{ \begin{aligned} \frac{1}{m'}, \quad i=l, \\ 0, \quad \quad i\neq l. \\ \end{aligned} \right. \end{eqnarray}$

(3.18)

From the orthogonality property of BPF, it is possible to expand functions into their Block Pulse series. This means that for every $f(t)\in{[0, 1)}$ we can write

$\begin{equation}\label{eq:29)} f(t) \approx \sum\limits_{i = 0}^{m'- 1} {{f_i}{b_i}(t) = {f^T}{B_{m'}}(t)}, \end{equation}$

(3.19)

where

$\begin{equation}\label{eq:30)} {f^T} = [{f_0}, {f_1}, \cdots, {f_{m'-1}}], \quad B_{m'}^T(t) = [{b_0}(t), {b_1}(t), \cdots, {b_{m'-1}}(t)], \end{equation}$

(3.20)

such that $f_i$ for $i=0, 1, 2, \cdots, m'-1$ are obtained by $f_i = m'\displaystyle\int_0^1 {f(t){b_i}(t)dt}$.

Similarly, the Chebyshev wavelets may be expanded into an ${m'}$-term block pulse functions (BPF) as

$\begin{equation}\label{eq:31)} \Psi (t) \approx {\Phi _{m' \times m'}}{B_{m'}}(t). \end{equation}$

(3.21)

We derive the Block Pulse operational matrix of the fractional integration $F^\alpha$ as following

$\begin{equation}\label{eq:32)} {I^\alpha }{B_{m'}}(t) = {F^\alpha }{B_{m'}}(t), \end{equation}$

(3.22)

where

$\begin{equation}\label{eq:33)} \begin{aligned} {F^\alpha } =\frac{1}{{{2^\alpha }{{m'}^\alpha }\Gamma (\alpha + 1)}}\left[{\begin{array}{ccccccccccccccc} 1&{{\xi _1}}&{{\xi _2}}&{{\xi _3}}& \cdots &{{\xi _{m'-1}}}\\ 0&1&{{\xi _1}}&{{\xi _2}}& \cdots &{{\xi _{m'-2}}}\\ 0&0&1&{{\xi _1}}& \cdots &{{\xi _{m'-3}}}\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots \\ 0&0&0&0& \cdots &{{\xi _1}}\\ 0&0&0&0&0&1 \end{array}} \right] \end{aligned} \end{equation}$

(3.23)

with ${\xi _k} = {(2k + 1)^\alpha } - {(2k - 1)^\alpha }, k=1, 2, \cdots, m'-1$.

Next, we derive the Chebyshev wavelet operational matrix of the fractional integration. Let

$\begin{equation}\label{eq:34)} {I^\alpha }\Psi (t) \approx {\rm P}_{m' \times m'}^\alpha \Psi (t). \end{equation}$

(3.24)

Using eqs. (3.21), (3.22}), we have

$\begin{equation}\label{eq:35)} {I^\alpha }\Psi (t) \approx {I^\alpha }{\Phi _{m' \times m'}}{B_{m'}}(t) = {\Phi _{m' \times m'}}{I^\alpha }{B_{m'}}(t) \approx {\Phi _{m' \times m'}}{F^\alpha }{B_{m'}}(t). \end{equation}$

(3.25)

From eqs. (3.24), (3.25), we get

$\begin{equation}\label{eq:36)} {\rm P}_{m' \times m'}^\alpha \Psi (t)\approx{\Phi _{m' \times m'}}{F^\alpha }{B_{m'}}(t). \end{equation}$

(3.26)

Then from eqs. (3.21), (3.26), the Chebyshev wavelet operational matrix of the fractional integration ${\rm P}_{m' \times m'}^\alpha$ is given by

$\begin{equation}\label{eq:37)} P_{m' \times m'}^\alpha = {\Phi _{m' \times m'}}{F^\alpha }\Phi _{m' \times m'}^{ - 1}. \end{equation}$

(3.27)

It should be noted that the operational matrix ${\rm P}_{m' \times m'}^\alpha$ contains many zero entries. This phenomenon makes calculations fast. The calculation for the matrix ${\rm P}_{m' \times m'}^\alpha$ is carried out once and is used to solve fractional order as well as integer order differential equations.

By approximating the function ${{D_ * }^vu(t)}$, we have

$\begin{equation}\label{eq:38)} {{D_ * }^vu(t)}u(t) \approx {C^T}\Psi (t) \end{equation}$

(4.1)

together with the initial states, we get

$\begin{eqnarray}\label{eq:39)} {D_ * }^{{\beta _i}}u(t) &\approx& {C^T}P_{m' \times m'}^{v - {\beta _i}}\Psi (t), i=1, 2, \cdots, v-1, \end{eqnarray}$

(4.2)

$\begin{eqnarray}\label{eq:40)} u(t) &\approx& {C^T}P_{^{m' \times m'}}^v\Psi (t) + \sum\limits_{k = 0}^{n - 1} {{u^{(k)}}(0)}\frac{t^k}{k!}. \end{eqnarray}$

(4.3)

Substituting eqs. (4.1), (4.2) and (4.3) into eq. (1.1), we have

$\begin{equation}\label{eq:41)} {C^T}\Psi (t) + \sum\limits_{i = 1}^{r - 1} {{\gamma _i}(t){C^T}} P_{m' \times m'}^{v - {\beta _i}}\Psi (t) + {\gamma _r}(t)({C^T}P_{m' \times m'}^v\Psi (t) + \sum\limits_{k = 0}^{n - 1} {{u^{(k)}}(0)\frac{t^k}{k!}) = g(t)}. \end{equation}$

(4.4)

Coefficient ${\gamma _i}(t) {(i=1, 2, \cdots, r)}$ can be dispersed into ${\gamma _i} (t_l)$ and $g(t)$ may be dispersed into $g(t_l) {(l=1, 2, \cdots, m')}$. Let

$\begin{eqnarray}\label{eq:42)} && {R_i} = \left[{\begin{array}{ccccccccccccccc} {{\gamma _i}({t_1})}&0& \cdots &0\\ 0&{{\gamma _i}({t_2})}& \cdots &0\\ \vdots&\vdots&\ddots&\vdots \\ 0&0& \cdots &{{\gamma _i}({t_{m'}})} \end{array}} \right] {(i=1, 2, \cdots, r)}, \end{eqnarray}$

(4.5)

$\begin{eqnarray} && \label{eq:43)} g = {\left[{\begin{array}{ccc} {g({t_1})-{\gamma _r}({t_1})\sum\limits_{k = 0}^{n-1} {{u^{(k)}}(0)\frac{t^k}{k!}} }& \cdots &{g({t_{m'}})-{\gamma _r}({t_{m'}})\sum\limits_{k = 0}^{n - 1} {{u^{(k)}}(0)\frac{t^k_{m'}}{k!}} } \end{array}} \right]^T}. \end{eqnarray}$

(4.6)

Therefore, eq. (4.4) can be written as

$\begin{equation}\label{eq:44)} \begin{aligned} {C^T}({\Phi _{m' \times m'}} + \sum\limits_{i = 1}^r {P_{m' \times m'}^{v - {\beta _i}}} {\Phi _{m' \times m'}}{R_i}) = {g^T}. \end{aligned} \end{equation}$

(4.7)

Eq. (4.7) is a linear system of algebraic equations.

Theorem 5.1  (see [18]) Let A function $f(t)$ define on the interval [0,1] may be expanded as

$\begin{equation}\label{eq:45)} f(t) = \sum\limits_{n = 1}^\infty {\sum\limits_{m = 0}^\infty {{c_{nm}}{\psi _{nm}}(t)} }, \tilde f(t) = \sum\limits_{n = 1}^{{2^{k - 1}}} {\sum\limits_{m = 0}^{M - 1} {{c_{nm}}{\psi _{nm}}(t)} }, \end{equation}$

(5.1)

then

$\begin{equation}\label{eq:46)} \sum\limits_{n = 1}^{{2^{k - 1}}} {\sum\limits_{m = 0}^{M - 1} {c_{nm}^2} } \le \int_0^1 {{f^2}(t)dt}, \end{equation}$

(5.2)

where

$\begin{equation}\label{eq:47)} {c_{nm}} = {(f(t), {\psi _{nm}}(t))_{{\omega _n}}} = \int_0^1 {{\omega _n}(t){\psi _{nm}}(t)f(t)dt}. \end{equation}$

(5.3)

Theorem 5.2  (see [18]) Let A function $f(t)$ be L²[0,1], and

$\begin{equation}\label{eq:48)} {R_{K, M}} = f(t) - \tilde f(t), \end{equation}$

(5.4)

then

$\begin{equation}\label{eq:49)} \lim\limits_{K, M \rightarrow \infty} \left\| {{R_{K, M}}} \right\| = 0, \end{equation}$

(5.5)

where $f(t), \tilde f(t)$ be defined as above and $K = {2^{k - 1}}$.

Also, we present an error estimation for eq. (5.4) with the residual error function. Let $u(t), \tilde{u(t)}$ be the exact solution and the Chebyshev wavelet solution of the problem (1.1) and (1.2), $R(t)=u(t)-\tilde{u(t)}$. Therefore, $R(t)$ satisfies the following problem

$\begin{eqnarray}\label{eq:50)} &&{D_ * }^vR(t) + \sum\limits_{i = 1}^{r - 1} {{\gamma _i}(t){D_ * }^{{\beta _i}}R(t)} + {\gamma _r}(t)R(t) = g(t) - h(t), \end{eqnarray}$

(5.6)

$\begin{eqnarray}\label{eq:51)} &&{R^{(i)}}(0) = 0, i = 0, 1, \cdots, n - 1, \end{eqnarray}$

(5.7)

where $h(t) = {C^T}\Psi (t) + \sum\limits_{i = 1}^{r - 1} {{\gamma _i}(t){C^T}P_{m' \times m'}^{v - {\beta _i}}\Psi (t)} + {\gamma _r}(t){C^T}P_{m' \times m'}^\alpha \Psi (t)$.

By solving the error problem (5.6)--(5.7) in the same way as Section 4, we get the approximation

$\begin{equation}\label{eq:52)} \tilde R(t) = \sum\limits_{n = 1}^{{2^{k - 1}}} {\sum\limits_{m = 0}^{M - 1} {c_{nm}^*{\psi _{nm}}(t)} } \end{equation}$

(5.8)

to $R(t)$. Here, the coefficients $c_{nm}^*, n = 1, 2, \cdots, {2^{k - 1}}, m = 0, 1, 2, \cdots, M - 1$ are determined by solving the error problem (5.6)--(5.7). Hence, the maximum absolute error can be estimated approximately by using $\max \left| {\tilde R(t)} \right|$. If the exact solution of the problem is not known, this error estimation can be used to test the reliability of the results.

Consider the equation

$\begin{equation}\label{eq:53)} {D^2}u(t) + \sin (t){D_{*}^{\frac{1}{2}}}u(t) + tu(t) = f(t), u(0) = u'(0) = 0, \end{equation}$

(6.1)

where $f(t) = {t^9} - {t^8} + 56{t^6} - 42{t^5} + \sin (t)(\frac{{32768}}{{6435}}{t^{\frac{{15}}{2}}} - \frac{{2048}}{{429}}{t^{\frac{{13}}{2}}})$.

One can easily check that $u(t)=t^8-t^7$ is the unique analytical solution. The comparison between approximate and exact solution for $k=3, M=5$ is presented in Fig. 1 and we list the absolution errors for ${M=3}$ and different values of $k$ and Ref. [Li and Zhao (2010)] (see[19]). From Table 1, we can achieve a better approximation by the Chebyshev wavelet method than Ref. [Li and Zhao (2010)]. We may also see that the error is smaller and smaller when $k$ increases. Therefore for better results, using a larger $k$ is recommended. The computational results show that the method in this article can be effectively used in numerical calculus for fractional differential equation with variable coefficient, and the method is also feasibility to the realistic fractional differential equation.

This article adopts the Chebyshev wavelet method to solve a class of multi-order fractional differential equations with variable coefficients. We derive the Chebyshev wavelet operational matrix of fractional order integration and use the wavelet basis together with operational matrix to reduce the factional differential equation to a system of algebraic equations. Furthermore, we present an error estimation for eq. (5.4) with the residual error function. The matrix elements of the discrete operators are provided explicitly and this in turn greatly simplifies the steps for obtaining solutions. An example is given to demonstrate that the method is effective and accurate for solving multi-order FDEs. It is obvious that the accuracy improves when we take a relatively small fixed ${M}$ and only increase ${k}$. Usually, it can reach the higher precision, even though ${M}, k$ are small.