\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small}
\vspace{9mm}}

\begin{document}
\title[\hfil fractional boundary value problems]
{Chebyshev finite difference method for fractional boundary value problems}

\author[H. Azizi, G. B. Loghmani \hfilneg]
{Hadi Azizi, Ghasem Barid Loghmani}  % in alphabetical order

\address{Hadi Azizi  \newline
Department of Mathematics, Taft Branch, Islamic Azad University, Taft, Iran.}
\email{hadiazizi1360@yahoo.com}

\address{Ghasem Barid Loghmani  \newline
Department of Mathematics, Yazd University, Yazd, Iran}
\email{loghmani@yazduni.ac.ir}

%\thanks{Submitted July 15, 2013. Published August 7, 2013.}
\subjclass[2000]{34A08}
\keywords{Chebyshev polynomials; Gauss-Lobatto points; Fractional differential  \hfill\break\indent equation; 
Finite difference}

\begin{abstract}
 This paper present a numerical method for fractional differential
equations using Chebyshev finite difference method. The fractional
derivatives are described in the Caputo sense. Numerical results
show that this method is of high accuracy and is more convenient and efficient for solving boundary value problems involving
fractional ordinary differential equations.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}\label{intro}

The idea of a derivative which interpolates between the
familiar integer order derivatives was introduced many years ago
and has gained increasing importance only in recent years due to
the development of mathematical models of a certain situations in
engineering, materials science, control theory, polymer modelling
etc. For example see \cite{magin,meral,old,old2,pod}.\\
\indent Most fractional order differential equations describing
real life situation, in general do not have exact analytical
solutions. Several numerical and approximate analytical methods
for ordinary differential equation have been extended to solve
fractional order differential equations. According to \cite{muj},
these methods include Adomian decomposition method \cite{ado},
homotopy perturbation method \cite{odi}, finite difference method
\cite{zhan}, fractional linear multi-step method \cite{lub},
extrapolation method \cite{kai} and predictor-corrector method
\cite{kai1}. Whilst much work has been published on fractional
differential equation to date, most of the work is focused on
initial value problems. Boundary value problems for fractional
differential equations has received attention only recently in
\cite{ahmad,bai,muj,shu}.\\
\indent In numerical analysis, Taylor polynomials are often used
to expand functions. Chebyshev polynomials can also be used to
expand functions. An advantage of using Chebyshev polynomials to
expand functions is the good representation of smooth function by
finite Chebyshev expansion if the function $y(t)$ is infinitely
differentiable \cite{elb}. Chebyshev finite difference method
(ChFD) has been successfully used in the numerical solution of
boundary value problems, boundary layer equations, nonlinear
system of second-order boundary value problems and Fredholm
integro-differential equations \cite{elb,elb1,saadat,dehghan}.\\
\indent Since the Chebyshev finite difference method for solving
ordinary differential equations is an appropriate method, we have
extended it to the fractional differential equations. So, the
objective of this work is to use ChFD method
to solve boundary value problems involving fractional differential equations.

\section{Preliminaries}

To develop our proposed technique, we require the
following definitions of
fractional integral, fractional derivatives and Chebyshev polynomials.\\
\textbf{Definition 2.1.}\cite{pod} Let $ \alpha \in R^{+} $. The
operator $ J_{a}^{\alpha} $, defined on the usual Lebesgue space $ L_{1}\left[ a,b \right]  $ by
$$
\nonumber J_{a}^{\alpha} f(x)=\dfrac{1}{\Gamma(\alpha)} \int
_{a}^{x}(x-t)^{\alpha -1} f(t) dt ,
$$
$$ \hspace*{-3.75cm} J_{a}^{0} f(x)= f(x) , $$
for $ a\leq x \leq b $, is defined to be the Riemann-Liouville
fractional
integral operator of order $ \alpha $ . Some properties of the operator $ J_{a}^{\alpha} $ are as follows:\\
 For $ f \in L_{1}\left[
a,b \right], \alpha , \beta \geq 0 $ and $\gamma >-1 $,
 \begin{flushleft}
$ (1)\quad ~~~ J_{a}^{\alpha}f(x) \hspace*{.1cm} exists
\hspace*{.1cm} for \hspace*{.1cm} almost \hspace*{.1cm} every
\hspace*{.1cm} x\in \left[ a,b \right] ,  $
\end{flushleft}
%
\begin{flushleft}
$ (2)\quad ~~~ J_{a}^{\alpha} J_{a}^{\beta} f(x)= J_{a}^{\alpha
+\beta} f(x), $
\end{flushleft}
%
\begin{flushleft}
$ (3)\quad ~~~ J_{a}^{\alpha} J_{a}^{\beta} f(x)= J_{a}^{\beta}
J_{a}^{\alpha} f(x), $
\end{flushleft}
%
\begin{flushleft}
$ (4)\quad ~~~ J_{a}^{\alpha} (x-a)^{\gamma} =
\dfrac{\Gamma(\gamma+1) }{\Gamma (\alpha +\gamma +1)}
(x-a)^{\alpha +\gamma}. $
\end{flushleft}
\textbf{Definition 2.2.}\cite{pod} The fractional derivative of $
f(x) $ in the Riemann-Liouville sense can be defined by
$$
D _{a}^{\alpha} f(x)=D^{m} J_{a}^{m-\alpha}f(x)=
\dfrac{d^{m}}{dx^{m}} \dfrac{1}{\Gamma (m-\alpha)} \int_{a} ^{x}
(x-t)^{m-\alpha -1} f(t) dt.
$$
Here $ m \in N $ satisfies the relation $ m-1 < \alpha \leq m
$, and $ f \in L_{1}\left[ a,b \right] $ .\\
For $ m-1 < \alpha ,\beta \leq m $ $, x>a $ and $ \gamma >-1 $ it
is know that:
\begin{flushleft}
$ (1)\quad ~~~ D_{a}^{\alpha} (x-a)^{\gamma} =\dfrac{\Gamma
(\gamma +1)}{\Gamma(\gamma-\alpha +1)} (x-a)^{\gamma -\alpha} , $
\end{flushleft}
\begin{flushleft}
$ (2)\quad ~~~ D_{a}^{\alpha} J_{a}^{\alpha} f(x)=f(x). $
\vspace*{.5cm}
\end{flushleft}
\textbf{Definition 2.3.}\cite{pod} The fractional derivative of $
f(t) $ in the Caputo sense can be defined by
$$
  { }_aD_{\ast}^{\alpha} f(t)=
\dfrac{1}{\Gamma (m-\alpha)} \int _{a}^{t} (t-s)^{m-\alpha -1}
f^{(m)} (s) ds.
$$
Here $ m-1 <\alpha \leq m , m \in N , t>a $ .\\
If $ m-1 < \alpha \leq m $, $ t>a $ then the following properties
hold:
\begin{flushleft}
$ (1)\quad ~~~ { }_aD_{\ast}^{\alpha} k=0, \vspace*{.25cm} $\\
$ (2)\quad ~~~ { }_aD_{\ast}^{\alpha} (J_{a}^{\alpha}f(t))=f(t) , \vspace*{.25cm} $\\
$ (3)\quad ~~~  J_{a}^{\alpha} ({
}_aD_{\ast}^{\alpha})=f(t)-\mathop{\sum}\limits _{k=0}^{m-1}
f^{(k)}(a) \dfrac{(t-a)^{k}}{k!}. $
\end{flushleft}
\textbf{Definition 2.4.}\cite{mason} The Chebyshev polynomials of
the first kind of degree $n$ are defined on the interval $[-1,1]$
by
$$
T_{n}(t)=\cos(n\arccos(t)),
$$
$T_{0}(t)=1$, $T_{1}(t)=t$ and they satisfy the recurrence
relations:
$$ T_{n+1}(t)=2tT_{n}(t)-T_{n-1}(t), \ \ n=1,2,\cdot\cdot\cdot.$$
\cite{clen} has given the following approximation of the function
$y(t)$
\begin{equation}\label{eq1}
 y(t)=\sum_{n=0}^{N}{^{''}r_{n}T_{n}(t)},
\end{equation}
where
$$r_{n}=\frac{2}{N}\sum\limits_{j = 0}^N {^{''}y(t_{j})T_{n}(t_{j})}, $$
and $ t_{j}=\cos(\frac{j\pi}{N}), \ \ j=0,1,...,N$ are
Gauss-Lobatto points. As in \cite{elb} the summation symbol with
double primes
denotes a sum with both the first and last terms halved.\\
The derivatives of the Chebyshev function are formed as follows:
\begin{equation}\label{eq2}
T_{n}'(t)=\sum_{\stackrel{k=0}{(n+k)odd}}^{n-1}\frac{2n}{c_{k}}T_{k}(t),
\end{equation}
\begin{equation}\label{eq3}
T_{n}''(t)=\sum_{\stackrel{k=0}{(n+k)even}}^{n-2}\frac{n}{c_{k}}(n^{2}-k^{2})T_{k}(t),
\end{equation}
\begin{equation}\label{eq4}
T_{n}'''(t)=\sum_{\stackrel{k=0}{(n+k)even}}^{n-2}\sum_{\stackrel{j=0}{(j+k)odd}}^{k-1}\frac{2kn}{c_{k}c_{j}}(n^{2}-k^{2})T_{j}(t),\\
\end{equation}
~~~~~~~~~~~~~~~~~~~~~~~~\vdots\\

where $c_{0}=2$ and $c_{i}=1$ for $i\geq1$.\\
\indent From equation (\ref{eq2})-(\ref{eq4}) it can be obtained
\begin{equation}\label{eq5}
y'(t)=\frac{4}{N}\sum_{n=0}^{N}{^{''}}\sum_{j=0}^{N}{^{''}\sum_{\stackrel{k=0}{(n+k)odd}}^{n-1}\frac{n}{c_{k}}y(t_{j})T_{n}(t_{j})T_{k}(t)},
\end{equation}
\begin{equation}\label{eq6}
y''(t)=\frac{2}{N}\sum_{n=0}^{N}{^{''}}\sum_{j=0}^{N}{^{''}\sum_{\stackrel{k=0}{(n+k)even}}^{n-2}\frac{n}{c_{k}}(n^{2}-k^{2})y(t_{j})T_{n}(t_{j})T_{k}(t)},
\end{equation}
\begin{equation}\label{eq7}
y'''(t)=\frac{4}{N}\sum_{n=0}^{N}{^{''}}\sum_{j=0}^{N}{^{''}\sum_{\stackrel{k=0}{(n+k)even}}^{n-2}\sum_{\stackrel{i=0}{(i+k)odd}}^{k-1}\frac{kn}{c_{i}c_{k}}(n^{2}-k^{2})y(t_{j})T_{n}(t_{j})T_{i}(t)}.
\end{equation}
~~~~~~\vdots\\

\indent The derivatives of the function $y(t)$ at the points
$t_{k}$ are given by \cite{dehghan,elb,ghaly}
\begin{equation}\label{eq8}
y^{(n)}(t_{k})=\sum_{j=0}^{N}d_{k,j}^{(n)} \ y(t_{j}), \ \ n=1,2,3,...
\end{equation}
where
\begin{equation}\label{eq9}
d_{k,j}^{(1)}=\frac{4\mu_{j}}{N}\sum_{n=0}^{N}\sum_{\stackrel{l=0}{(n+l)odd}}^{n-1}\frac{n\mu_{n}}{c_{l}}T_{n}(t_{j})T_{l}(t_{k}),
\ \ k,j=0,1,...,N,
 \end{equation}
\begin{equation}\label{eq10}
d_{k,j}^{(2)}=\frac{2\mu_{j}}{N}\sum_{n=0}^{N}\sum_{\stackrel{l=0}{(n+l)even}}^{n-2}\frac{n\mu_{n}(n^{2}-l^{2})}{c_{l}}T_{n}(t_{j})T_{l}(t_{k}),
\ \ k,j=0,1,...,N,
 \end{equation}
 \begin{equation}\label{eq11}
d_{k,j}^{(3)}=\frac{4\mu_{j}}{N}\sum_{n=0}^{N}\sum_{\stackrel{l=0}{(n+l)even}}^{n-2}\sum_{\stackrel{i=0}{(i+l)odd}}^{l-1}\frac{nl\mu_{n}(n^{2}-l^{2})}{c_{i}c_{l}}T_{n}(t_{j})T_{i}(t_{k}),
\ \ k,j=0,1,...,N,
 \end{equation}
 ~~~~\vdots\\

\hspace{-0.6 cm}with $\mu_{0}=\mu_{N}=\frac{1}{2}$ and $\mu_{j}=1$ for
$j=1,2,...,N-1$.\\
 \indent As can be seen from equation (\ref{eq8}), the derivatives of $y(t)$ at any point from the
 Gauss-Lobatto nodes are expanded as a linear combination of the values
 of the function at these points.

\section{Chebyshev finite difference method}
\indent In this section by using ChFD method the boundary value
problem of fractional order of the form
\begin{equation}\label{eq12}
{
}_{-1}D_{\ast}^{\alpha}y(t)+a_{2}(t)y''(t)+a_{1}(t)y'(t)+a_{0}(t)y(t)=f(t),
\ \ \ -1\leq t\leq 1
\end{equation}
with the boundary condition
\begin{equation}\label{eq13}
y(-1)=\beta_{1}, \ \ y(1)=\beta_{2}
\end{equation}
where $1<\alpha\le 2$, $a_{0}(t),a_{1}(t),a_{2}(t)$ and $f(t)$
are given function in $L^{2}[-1,1]$ and $\beta_{1},\beta_{2}$ are
real given numbers is solved.\\
 \indent In order to find the solution $y(t)$ in equation (\ref{eq12}), one needs to first
 calculate equation (\ref{eq12}) in Gauss-Lobatto nodes $t_{k}$ for
 $k=1,2,...,N-1$ and use equations (\ref{eq8}), (\ref{eq9}), (\ref{eq10}) as well as definition 2.3 to obtain
\begin{equation}\nonumber
\nonumber
\frac{1}{\Gamma(2-\alpha)}\int_{-1}^{t_{k}}(t_{k}-s)^{-\alpha}\sum_{p=0}^{N}{^{''}}\sum_{q=0}^{N}{^{''}
\sum_{\stackrel{r=0}{(p+r)even}}^{p-2}\frac{2p}{Nc_{r}}(p^{2}-r^{2})y(t_{q})T_{p}(t_{q})T_{r}(s)}ds\\
\end{equation}
\begin{equation}\label{eq14}
+\sum_{j=0}^{N}y(t_{j})a_{2}(t_{k})d_{k,j}^{(2)}+\sum_{j=0}^{N}y(t_{j})a_{1}(t_{k})d_{k,j}^{(1)}+a_{0}(t_{k})y(t_{k})=f(t_{k}).
\end{equation}
 For $k=0$ and $k=N$ the boundary condition (\ref{eq13}) is
used to obtain
\begin{equation}\label{eq15}
y(t_{0})=y(1)=\beta_{2}, \ \ y(t_{N})=y(-1)=\beta_{1}.
\end{equation}
Therefore equations (\ref{eq14}) and (\ref{eq15}) generate a set of $N+1$
algebraic equations, which can be solved for unknown
$y(t_{0}),y(t_{1}),...,y(t_{N-1})$ and $y(t_{N})$ and thus
$y(t)$ in equation (\ref{eq1}) can be calculated.\\

\hspace*{-0.4cm}\textbf{Remark 3.1.} For system of fractional differential equation we can apply Chebyshev finite difference method similarly. Therefore, we solve 
two examples in the next section.\\ 
\textbf{Remark 3.2.} If the interval under consideration in the
boundary condition is $[a,b]$, we transform it to the interval
$[-1,1]$ by the change of variable,
$$t=\frac{b-a}{2}x+\frac{b+a}{2}.$$ 

\section{Numerical examples} 
\indent In this section, we apply the ChFD method to solve different types of linear fractional differential equations to show the
accuracy of the proposed method. All results were computed using Maple 13. \\

4.1. Ordinary fractional differential equations\\

\textbf{Example 4.1.1.} Consider the boundary value problem of
fractional order 
$$
{
}_{-1}D_{\ast}^{1.4}y(t)+ty(t)=\frac{(t+1)^{\frac{3}{5}}}{4\Gamma(0.6)}(25t-15)+t^{4},
\ \ \ -1\leq t\leq1,
$$
$$
y(-1)=-1, \ \ \ y(1)=1.
$$
The exact solution of this problem is $y(t)=t^{3}$. We apply the
proposed method for this problem and achieve the
exact solution.\\

\textbf{Example 4.1.2.} \cite{muj} Consider the boundary value problem 
$$ 
{
}_0D_{\ast}^{\alpha}y(t)-\frac{1-t}{(1+t)^{2}}y(t)=\frac{1}{(1+t)^{2}}, 
\ \ 1<\alpha \le 2
$$ 
$$ 
y(0)=1, \ \ y(1)=\frac{1}{2}. 
$$
We apply the ChFD method and the approximate solution of this problem for different values of $\alpha$ are shown in figure  1.\\
When $\alpha=2$, the exact solution is $y(t)=\frac{1}{1+t}$ 
\cite{li}. For this case we compare our results with methods based 
on Taylor's expansion \cite{li} and Haar wavelet method 
\cite{muj}. The evaluated absolute errors between the exact
solution and numerical solution are shown in Table 1. 
\begin{figure}[h!] 
\begin{tabular}{c}
\includegraphics[width=7cm,height=5cm]{image/fig2.eps} 
\end{tabular}
\caption{ Approximate solution with N=8 and different values of $\alpha$ for example 4.1.2.  } 
\centerline{  \begin{footnotesize} \end{footnotesize}  }
\end{figure}
\begin{center} 
{\footnotesize Table 1:Comparison of the ChFD method for N=20 with
the Haar wavelet \cite{muj} and methods in ref. \cite{li}.}
{\small
\[\begin{array}{|c|cc|c|c|}\hline 
&~~~~~~~~~~~~~~~$Methods in ref. \cite{li}$& &$Haar wavlet \cite{muj}$&$ChFD method$ \\ \hline
t&|\tilde{y}(t)-y(t)|&|\hat{y}(t)-y(t)|&|y_{Haar}-y(t)|&|y_{CHFD}(t)-y(t)| \\\hline

0.1&1.590477\times10^{-4}&2.281673\times10^{-4}&1.855711\times10^{-7}&5.196667\times10^{-17}\\

0.2&1.047506\times10^{-4}&2.469834\times10^{-4}&1.307581\times10^{-7}&8.692174\times10^{-17} \\

0.3&2.784706\times10^{-5}&1.858469\times10^{-4}&5.238972\times10^{-8}&2.484844\times10^{-16} \\

0.4&4.550068\times10^{-5}&1.157510\times10^{-4}&2.279295\times10^{-8}&1.933657\times10^{-16} \\

0.5&1.167467\times10^{-4}&6.412753\times10^{-5}&6.420743\times10^{-9}&1.365740\times10^{-17}\\

0.6&1.812416\times10^{-4}&3.403860\times10^{-5}&2.6709901\times10^{-9}&1.654892\times10^{-16} \\

0.7&2.228086\times10^{-4}&2.153643\times10^{-5}&5.052727\times10^{-9}&1.822994\times10^{-16} \\

0.8&2.213083\times10^{-4}& 1.719699\times10^{-5}&1.285991\times10^{-8}&4.105945\times10^{-17}\\
0.9&1.543950\times10^{-4}&1.205234\times10^{-5}&2.474367\times10^{-8}&3.606577\times10^{-17} \\
\hline
\end{array}\]}
\end{center}
4.2. Multi-order fractional differential equations\\

\textbf{Example 4.2.1.} Consider the boundary value problem of
fractional order
$$
{
}_0D_{\ast}^{1.25}y(t)+y''(t)+y'(t)+ty(t)=\frac{32}{12\Gamma(0.75)}t^{0.75}+t^{3}+2t+2,
\ \ \ 0\leq t\leq1,
$$
$$
y(0)=0, \ \ y(1)=1,
$$
with the exact solution $y(t)=t^{2}$.\\
By using Chebyshev finite difference method with $N=2$, we obtain
the exact solution $y(t)=t^{2}$.\\

\textbf{Example 4.2.2.} Consider the boundary value problem of
fractional order
$$
{
}_0D_{\ast}^{1.5}y(t)+y'(t)+\frac{y(t)}{e^{t}}=e^{t}+erf(\sqrt{t})e^{t}+1,
\ \ \ 0\leq t \leq1,
$$
$$
y(0)=1, \ \ \ y(1)=e,
$$
with the exact solution $y(t)=e^{t}$.\\
Now we define the maximum errors for $y_N
(t)=\sum_{n=0}^{N}{^{''}r_{n}T_{n}(t)}$ as
$$
e_{N}=||y_{N}(t)-y_{exact}(t)||_{\infty}=max\{|y_{N}(t)-y_{exact}(t)|,\
0\leq t\leq1 \}
$$
where $y_{N}$ is the computed result with $N$ and $y_{exact}$ is
the exact solution. In Table 2 we give the errors $e_{N}$ for
different values of $N$.
\begin{center}
{\footnotesize Table 2: The maximum errors of $e_{N}$ for
different values of N for example 4.2.2}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
    %after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  \ \  N ~~~~& 4 & 6 & 8 & 10 & 12 \\
\hline
  \ \  $e_N ~~~~$ & $10^{ - 4}$ & $10^{ - 7}$ & $10^{ - 10}$ & $5\times 10^{ - 14}$ & $2.5\times 10^{ - 16}$ \\
  \hline
\end{tabular}
\end{center}
\vspace{0.5cm}
\textbf{Example 4.2.3.} \cite{muj} Consider the boundary value problem
$$
{ }_0D_{\ast}^{\alpha}y(t)={ }_0D_{\ast}^{\beta}y(t)-e^{t-1}-1, \
\ 1<\alpha \le2, \ \ 0<\beta \le1,
$$
$$
y(0)=0, \ \ y(1)=0.
$$
For the general case the exact solution of the problem is not
known. However, for $\alpha=2$ and $\beta=1$ the problem has the
exact solution $y(t)=t(1-e^{t-1}).$ This problem was solved
numerically in \cite{wang} using combined homotopy perturbation
method and Green function method. Also in \cite{muj} it was solved
by Haar wavelet method. We compare ChFD method and the two
mentioned methods in Table 3. Computer plots for $\beta=1$ and
different values of $\alpha$ given in Figure 2 show that as
$\alpha$ approaches to 2, the
corresponding solutions of fractional order differential equation approach to the solutions of integer order differential equation.\\

\begin{figure}[h!]
\begin{tabular}{c}
\includegraphics[width=7cm,height=5cm]{image/fig3.eps}
\end{tabular}
\caption{ {\small Approximate solution with N=8, $\beta=1$ and different values of $\alpha$ for example 4.2.3. }}
\centerline{  \begin{footnotesize} \end{footnotesize}  }
\end{figure}
\begin{center}
{\footnotesize {\small Table 3: Comparison of the ChFD method for N=20 with the Haar wavelet \cite{muj} and HPM \cite{wang}.}}
{\small
\[\begin{array}{ccccc}\hline
t~~~~~~&~~~~~~$HPM \cite{wang}$&~~~~~~$Haar wavlet \cite{muj}$&~~~~~~$exact solution$&~~~~~~|y_{CHFD}-y(t)|\\\hline
0.1~~~~~~&~~~~~~0.05934820&~~~~~~0.05934300&~~~~~~0.05934303&~~~~~~1.0213\times10^{-26}\\
0.2~~~~~~&~~~~~~0.11014318&~~~~~~0.11013418&~~~~~~0.11013421&~~~~~~1.0682\times10^{-26}\\
0.3~~~~~~&~~~~~~0.15103441&~~~~~~0.15102438&~~~~~~0.15102441&~~~~~~1.1182\times10^{-26}\\
0.4~~~~~~&~~~~~~0.18048329&~~~~~~0.18047531&~~~~~~0.18047535&~~~~~~1.1682\times10^{-26}\\
0.5~~~~~~&~~~~~~0.19673826&~~~~~~0.19673463&~~~~~~0.19673467&~~~~~~1.2201\times10^{-26}\\
0.6~~~~~~&~~~~~~0.19780653&~~~~~~0.19780792&~~~~~~0.19780797&~~~~~~1.2733\times10^{-26}\\
0.7~~~~~~&~~~~~~0.18142196&~~~~~~0.18142718&~~~~~~0.18142725&~~~~~~1.3297\times10^{-26}\\
0.8~~~~~~&~~~~~~0.14500893&~~~~~~0.14501532&~~~~~~0.14501540&~~~~~~1.3946\times10^{-26}\\
0.9~~~~~~&~~~~~~0.08564186&~~~~~~0.08564623&~~~~~~0.08564632&~~~~~~1.4753\times10^{-26}\\ \hline
\end{array}\]}
\end{center}
\vspace{0.5cm} \textbf{Example 4.2.4.} The Bagley-Torvik equation
is of the form
$$
ay''(t)+b~{ }_0D_{\ast}^{\alpha}y(t)+cy(t)=g(t) ,\ \ \  t\in[0,1]
$$
$$y(0)=0, \ \ y(1)=y_{0}$$
where $a,b,c\in \mathbb{R}$ and $a\neq0$. The Bagley-Torvik
equation which involves fractional derivative of order
$\frac{1}{2}$ or $\frac{3}{2}$ arises in the modelling of the
motion of a rigid plate in a Newtonian fluid and a gas in a fluid.
According to \cite{muj} several authors, such as Podlubny
\cite{pod}, Saha Ray and Bera \cite{saha}, Diethelm and Ford
\cite{kai2}, Diethelm et al. \cite{kai1} and Wang and Wang
\cite{zwang}, have proposed different techniques for the solution
of Bagley-Torvik equation. In
the present work, we solve this equation with two-point boundary conditions by using the Chebyshev finite difference method.\\
We choose $\alpha=\frac{3}{2}$, $a=b=c=1$, $y_{0}=1$ and
$g(t)=4\sqrt{\frac{t}{\pi}}+t^{2}+2$. It can be easily verified
that the exact solution is
$y(t)=t^{2}$. For solving this special case we applied the ChFD method with $N=2$ and obtained the exact solution.\\

4.3. System of fractional differential equations\\

\textbf{Example 4.3.1.}  Consider the system of fractional differential equations
 \[
\left\{ {{\begin{array}{*{20}c}
 {
}_0D_{\ast}^{1.5}y_{1}(t)+y'_{2}(t)+t^{2}y'_{1}(t)=f_{1}(t), \ \ 0\leq t\leq1 \hfill \\
 {
}_0D_{\ast}^{1.5}y_{2}(t)+(t-1)y''_{1}(t)=f_{2}(t),~~~ \ \ 0\leq t\leq1 \hfill \\
\end{array} }} \right.
\]
with the boundary conditions
$$
y_{1}(0)=0, y_{1}(1)=1, \ \ y_{2}(0)=y_{2}(1)=0,
$$
where $f_{1}(t)=\frac{4}{\Gamma(0.5)}\sqrt{t}+2t^{3}+3t^{2}-1$ and
$f_{2}(t)=\frac{8}{\Gamma(0.5)}t^{\frac{3}{2}}+2(t-1)$. The exact
solutions are $y_{1}(t)=t^{2}$ and $y_{2}(t)=t^{3}-t$. For solving
the above equations we put
$$
y_{i}(t)=\sum_{n=0}^{N}{^{''}r_{n}^{i}T_{n}(t)}, \ \ r_{n}^{i}=\frac{2}{N}\sum\limits_{j = 0}^N {^{''}y_{i}(t_{j})T_{n}(t_{j})}, \ i=1,2,
$$
by using ChFD method and boundary conditions we obtain $2N+2$
algebraic equations. By solving this system we obtain unknown
coefficients $y_{1}(t_{k})$ and $y_{2}(t_{k})$ for $k=0,1,...,N$.
Now if we put $N=3$ and using the ChFD method for solving the
above system of fractional differential equation we obtain the
exact solutions $y_{1}(t)=t^{2}$
and $y_{2}(t)=t^{3}-t$.\\

\textbf{Example 4.3.2} \cite{lu} Consider the system of fractional
differential equation
\[
\left\{ {{\begin{array}{*{20}c}
 {
}_0D_{\ast}^{\alpha}y_{1}(t)+(2t-1)y'_{1}(t)+cos(\pi t)y'_{2}(t)=f_{1}(t), \ \ 0\leq t\leq1 \hfill \\
 {
}_0D_{\ast}^{\beta}y_{2}(t)+ty_{1}(t)=f_{2}(t),~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \ \ 0\leq t\leq1 \hfill \\
\end{array} }} \right.
\]
with the boundary conditions
$$
y_{1}(0)=0, y_{1}(1)=0, \ \ y_{2}(0)=y_{2}(1)=0,
$$
where $1<\alpha, \beta \leq 2$ and $f_{1}(t), f_{2}(t)$ are
determined by substituting exact solutions in the above equations.
For $\alpha=\beta=2$ the system have exact solutions
$y_{1}(t)=sin(\pi t)$ and $y_{2}(t)=t^{2}-t$. This problem was
solved in \cite{lu} using variational iteration method. We compare
the absolute error for $y_{1}(t)$ in Table 4. Note that by using
VIM and ChFD method we can obtain the exact solution for
$y_{2}(t)$. Now we suppose $\alpha=1.5, \beta=1.2$ and the exact
solutions for this case is the same as above. We apply ChFD method
and obtain the exact solution for $y_{2}(t)$ and for $y_{1}(t)$
the absolute errors are reported in Table 5.
\begin{center}
\hspace{-0.9cm} {\footnotesize  Table 4: Comparison of the CHFD method for N=16 with the VIM \cite{lu}
.}
{\small
\[\begin{array}{ccc}\hline
t&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $VIM \cite{lu}$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &$CHFD method$\\\hline
0.1&0.0003&8.1492\times10^{-17}\\
0.2&0.0025&6.5735\times10^{-18}\\
0.3&0.0078&5.0842\times10^{-18}\\
0.4&0.0166&1.2769\times10^{-18}\\
0.5&0.0277&3.5625\times10^{-18}\\
0.6&0.0387&1.2770\times10^{-18}\\
0.7&0.0459&5.0842\times10^{-18}\\
0.8&0.0449&6.5832\times10^{-18}\\
0.9&0.0309&8.1498\times10^{-17}\\ \hline
\end{array}\]}
\end{center}
\vspace{0.5cm}
\begin{center}
\hspace{-0.9cm} {\footnotesize  Table 5: Absolute error for $\alpha=1.5, \beta=1.2$ with N=16.}
{\small
\[\begin{array}{cc}\hline
t\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &|y_{CHFD}-y_{1}(t)|\\\hline 
0.1\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &4.3004\times10^{-16}\\ 
0.2\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &7.9298\times10^{-16}\\ 
0.3\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &1.1305\times10^{-15}\\ 
0.4\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &1.2186\times10^{-15}\\ 
0.5\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &1.2259\times10^{-15}\\ 
0.6\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &1.2196\times10^{-15}\\ 
0.7\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &8.9839\times10^{-16}\\ 
0.8\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &6.6944\times10^{-16}\\ 
0.9\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ &2.1484\times10^{-16}\\ \hline 
\end{array}\]}
\end{center}

\section{Conclusion}
\indent In this paper we have developed a numerical scheme for
solving fractional boundary value problems. Our proposed scheme
was based on the Chebyshev finite difference method. The results
obtained were compared with exact solutions and the results agree
well with exact solutions even for small values of $N$. In some
examples, the method achieved the exact solutions. The results
obtained indicate that this approach can solve boundary value
problems involving fractional ordinary differential equation
effectively.

\begin{thebibliography}{}

\bibitem{ahmad} B. Ahmad, Juan J. Nieto, \emph{Existence of solutions for nonlocal boundary value problems of higher-order nonlinear fractional differential equations
}, Abstr.
Appl. Anal., vol. 2009, Article ID 494720, pages 9 $<doi:10.1155/2009/494720>$.

\bibitem{bai} Z. Bai, \emph{On positive solutions of a nonlocal fractional boundary value problem}, Nonlinear Anal., \textbf{72}, 916-924, (2010)

\bibitem{clen} C.W. Clenshaw, A.R. Curtis, \emph{A method for numerical
integration on an automic computer}, Numer. Math., \textbf{2}, 197-205, (1960).

\bibitem{ado} V. Daftardar-Gejji, H. Jafari, \emph{Solving a multi-order fractional differential equation using adomian decomposition}, Appl. Math. Comput., \textbf{189}, 541-548,(2007).

\bibitem{dehghan} M. Dehghan, A. Saadatmandi, \emph{Chebyshev finite difference method for Fredholm
integro-differential equation}, International Journal of Computer Mathematics, \textbf{85}, 123-130, (2008).

\bibitem{kai2} K. Diethelm, N.J. Ford, \emph{Numerical solution of the Bagley-Torvik equation}, BIT \textbf{42}, 490-507, (2002).

 \bibitem{kai1} K. Diethelm, N.J. Ford, A.D. Freed, \emph{A predictor-corrector approach for the numerical solution of fractional differential equations}, Nonlinear Dyn. \textbf{29}, 3-22, (2002).


\bibitem{kai} K. Diethem, Guido walz, \emph{Numerical solution of fractional order differential equation by extrapolation}, Numerical Algorithms, 231-253 (1997).

\bibitem{elb} E.M.M Elbarbary, \emph{Chebyshev finite difference approximation for the boundary value problems}, Appl. Math. Comput. \textbf{139}, 513-523, (2003).

\bibitem{elb1} E.M.M Elbarbary, \emph{Chebyshev finite difference method for the solution of boundary-layer equation}, Appl. Math. Comput., \textbf{160}, 487-498, (2005).

\bibitem{ghaly} A.Y. Ghaly, M.A. Seddeek, \emph{Chebyshev finite difference method for the effects of chemical reaction, heat and mass transfer on laminar
flow along a semi infinite horizontal plate with temperature dependent viscosity}, Chaos, Solitons and Fractals, \textbf{19}, 61-70, (2004).

\bibitem{lu} Lu Junfeng. \emph{Variational iteration method for solving a nonlinear system of second-order boundary value problems}, Comput. Math. Appl., \textbf{54}, 1133-1138, (2007).

\bibitem{lub} C. Lubich, \emph{Fractional linear multistep methods for Abel-Volterra integral equations of the second kind}, Math. Comput., \textbf{45}, 463-469, (1985).

\bibitem{magin} R.L. Magin, \emph{Fractional calculus models of complex dynamics in biological tissues}, Comput. Math. Appl., \textbf{59}, 1586-1593, (2010).

\bibitem{mason} J.C. Mason, D.c. \emph{Handscomb, Chebyshev
polynomials}, CRC press, Boca Raton, (2003).

\bibitem{meral} F.C. Meral, T.J. Royston, R. Magin, \emph{Fractional calculus in viscoelasticity: an experimental study}, Commun. Nonl. Sci. Num. Sim., \textbf{15}, 939-945, (2010).

\bibitem{muj} Mujeeb ur Rehman, Rahmat Ali Khan, \emph{A numerical method for solving boundary value problems for
fractional differential equations}, Appl. Math. Modell., \textbf{36}, 894-907, (2012).

\bibitem{odi} Z. Odibat, S. Momani, V. Suat Erturk, \emph{Generalized differential transform method: application to differential equations of fractional order}, Appl. Math. Comput., \textbf{197}, 467-477, (2008).

\bibitem{old} K.B. Oldham, J.Spanier, \emph{The fractional calculus}, Academic press, San Diego, (1999).

\bibitem{old2} K.B. Oldham, \emph{Fractional differential equations in electrochemistry}, Adv. Eng. Soft., \textbf{41}, 9-12, (2010).

\bibitem{pod} I. Podlubny, \emph{Fractional differential equation}, Academic press, New
yourk, (1999).

\bibitem{saadat} A. Saadatmandi, J. A. Farsangi, \emph{Chebyshev finite difference method for a nonlinear system of second-order boundary value problems}, Appl. Math. Comput., \textbf{192}, 586-591, (2007).

\bibitem{saha} S. Saha Ray, R.K. Bera, \emph{Analytical solution of the Bagley-Torvik equation by Adomian decomposition method}, Appl. Math. Comput., \textbf{168}, 398-410, (2005).

\bibitem{shu} Z. Shuqin, \emph{Existence of solution for boundary value problem of fractional order}, Acta Math. Sci., \textbf{26}, 220-228, (2006).

\bibitem{wang} Y. Wang, H. Song, D. Li, \emph{Solving two-point boundary value problems using combined homotopy perturbation method and Greens function method},
Appl. Math. Comput., \textbf{212}, 366-376, (2009).

\bibitem{zwang} Z.H. Wang, X. Wang, \emph{General solution of the Bagley-Torvik equation with fractional-order derivative}, Commun. Nonlinear. Sci. Numer. Simul., \textbf{15}, 1279-1285, (2010).

\bibitem{li} Xian-Fang Li, \emph{Approximate solution of linear ordinary differential equations with variable coefficients}, Math. Comput. Simul., \textbf{75}, 113-125, (2007).

\bibitem{zhan} Y. Zhang, \emph{A finite difference method for fractional partial differential equation}, Appl. Math. Comput., \textbf{215}, 524-529, (2009).

\end{thebibliography}

\end{document}