\documentclass[11pt,twoside]{article}
\usepackage{amsmath, amsthm, amscd, amsfonts, amssymb, graphicx, color}
\usepackage[bookmarksnumbered, colorlinks]{hyperref} \usepackage{float}
\usepackage{lipsum}
\usepackage{afterpage}
\usepackage[labelfont=bf]{caption}
\usepackage[nottoc,notlof,notlot]{tocbibind} 
%\renewcommand\bibname{References}
\def\bibname{\Large \bf  References}
\usepackage{lipsum}
\usepackage{fancyhdr}
\pagestyle{fancy}
\fancyhf{}
\renewcommand{\headrulewidth}{0pt}
\fancyhead[LE,RO]{\thepage}
\thispagestyle{empty}
%\afterpage{\lhead{new value}}

\fancyhead[CE]{F. AUTHOR, S. AUTHOR AND T. AUTHOR} 
\fancyhead[CO]{ENTER THE TITLE OF THE PAPER ALL CAPITAL}



%\topmargin=-1.6cm
\textheight 17.5cm%
\textwidth  12cm %
\topmargin   8mm  %
\oddsidemargin   20mm   %
\evensidemargin   20mm   %
\footskip=24pt     %

\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{xca}[theorem]{Exercise}
%\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\renewenvironment{proof}{{\bfseries \noindent Proof.}}{~~~~$\square$}
\makeatletter
\def\th@newremark{\th@remark\thm@headfont{\bfseries}}
\makeatletter



  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If you want to insert other packages. Insert them here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\long\def\symbolfootnote[#1]#2{\begingroup%
%\def\thefootnote{\fnsymbol{footnote}}\footnote[#1]{#2}\endgroup}



 \def \thesection{\arabic{section}}
 

\begin{document}
%\baselineskip 9mm
%\setcounter{page}{}
\thispagestyle{plain}
{\noindent Journal of Mathematical Extension \\
Journal Pre-proof}\\
ISSN: 1735-8299\\
URL: http://www.ijmex.com\\
Original Research Paper\\
\vspace*{9mm}

\begin{center}

{\Large \bf 
A finite difference approximation for the solution of the space fractional diffusion equation\\}
%{\bf Do You Have a Subtitle? \\ If so, Write It Here} 


\let\thefootnote\relax\footnote{\scriptsize Received: January 2008; Accepted: February 2009 }


{\bf  H. Hajinezhad$^*$\let\thefootnote\relax\footnote{$^*$Corresponding Author}}\vspace*{-2mm}\\
\vspace{2mm} {\small  Payame Noor Univesity} \vspace{2mm}

{\bf A. R. Soheili}\vspace*{-2mm}\\
\vspace{2mm} {\small  Ferdowsi University of Mashhad} \vspace{2mm}
\end{center}

\vspace{4mm}


{\footnotesize
\begin{quotation}
{\noindent \bf Abstract.} The objective of this paper is to present a finite difference scheme that estimates the space fractional diffusion equation with the Caputo fractional derivative. The proposed scheme's stability, and convergence are proved. In order to evaluate the effectiveness of this scheme, a series of tests are conducted. The results of these tests demonstrate the reliability and accuracy of the proposed scheme.
\end{quotation}
\begin{quotation}
\noindent{\bf AMS Subject Classification:} 65M06; 65M12 ; 26A33; 35R11 

\noindent{\bf Keywords and Phrases:} Space fractional diffusion equation, Convergence, Stability, Finite difference method
\end{quotation}}

\section{Introduction}
In this paper, the space fractional diffusion equation
\begin{equation} \label{1}
\begin{aligned} 
\frac{\partial u(x,t)}{\partial t}=\alpha(x)\frac{\partial^\beta u(x,t)}{\partial x^ \beta}+f(x,t), \ 0<x<L,\ 0<t\leqslant T, \ 1< \beta <2
\end{aligned}
\end{equation} \\
with initial condition
\begin{align}
&u(x,0)=\Psi_{t_{0}} (x), \quad \quad 0\leqslant x\leqslant L, \label{2}
\end{align}
and boundary conditions
\begin{align}
&\frac{\partial u(0,t)}{\partial x}=\Psi_{x_{0}}(t), \quad 0 \leqslant t \leqslant T, \label{3}\\
&u(L,t)=\Psi_{x_{1}}(t), \quad \quad 0\leqslant t \leqslant T, \label{4}
\end{align}
is approximated, where $ u(x,t) $ is an unknown function and the space fractional derivative is assumed to be based on the Caputo fractional derivative as follows.
\begin{align}
\frac{\partial^\beta u(x,t)}{\partial x^ \beta}=\frac{1}{\Gamma(2-\beta)} \int_{0}^{x}{\frac{\partial^2 u(s,t)}{\partial s^2} {(x-s)}^{1-\beta} ds}, \quad 1<\beta <2. \label{5}
\end{align}\\

The space fractional diffusion equation with the Riemann fractional derivative was approximated by a shifted $Gr\ddot{u}nwald$ finite difference formula in \cite{meerschaert2006}. In \cite{tadjeran2006}, for equation \eqref{1} with the Riemann fractional derivative, the Crank--Nicolson method was applied based on a $Gr\ddot{u}nwald$ formula then an extrapolation was used to obtain a second-order approximation. \\ 

Equation \eqref{1} with fractional derivative \eqref{5} was solved numerically using orthogonal polynomials by some authors. The Legendre polynomials with the tau method were used to approximate this equation in \cite{saadatmandi2011}. Ren et al. \cite{ren2013} applied the shifted Chebyshev polynomials with the tau method to obtain an approximation for this equation. Some other authors used Chebyshev polynomials to estimate the solution of equation \eqref{1} with fractional derivative \eqref{5} \cite{khader2011}, \cite{safdari2019}, \cite{kheybari2019}, \cite{agarwal2018}. Khader \cite{khader2011} applied the Chebyshev polynomials to reduce this equation to the system of ordinary differential equations and then a finite difference approximation was used to obtain the numerical solution of this system. Safdari et al. \cite{safdari2019} approximated this equation by using the compact finite difference to obtain the semi-discretization in the time derivative and then used the Chebyshev collocation method to estimate the space fractional derivative. \\

In this paper, we propose a novel finite difference method for approximating equation \eqref{1} in Caputo sense \eqref{5} subject to conditions \eqref{2}-\eqref{4}. Our method has a distinct advantage over other methods used for space fractional diffusion with the Caputo derivative. Specifically, we demonstrate that our proposed scheme is unconditionally stable and convergent through rigorous proof. To evaluate the accuracy of our method, we conduct several numerical tests.\\

The structure of our paper is as follows:  The discretization of equation \eqref{1} is explained in the next Section. Section \ref{s-c} is dedicated to proving the stability and convergence of our proposed scheme. In Section \ref{test}, the numerical tests are provided. The final Section presents the conclusion.

    
\section{Finite difference method for the problem } \label{discrete}
The discretization of equation \eqref{1} using a proposed finite difference method is explained in this section.\\

For the finite difference scheme, assume the grid sizes in time and space as $ \Delta t $ and $ \Delta x $, respectively. Then, $ x_{j}=j \Delta x \ (j=0, 1, ..., J)$ and $ t^{n}=n \Delta t \ (n=0, 1, ..., N)$, where
$ J \Delta x=L $ and $ N \Delta t=T $. Assume $ u_{j}^{n} $ is the value of $ u(x_{j},t^{n}) $ for $ j=0, 1, ..., J$ and $\ n=0, 1, ..., N $. \\

The following lemma provides the essential tools for the discretization of equation \eqref{1}.
\begin{lemma} \label{lemma1}
Assume $ a_{s}=(s+1)^{2-\beta}-s^{2-\beta}, \ (s=0, 1, ..., \ 1< \beta <2) $. Then, the discretization of $ \frac{\partial ^{\beta} u(x,t)}{\partial x ^\beta} $ at $ (x_{j}, t^{n}) $ for $1 \leqslant j \leqslant J-1 $ and $0\leqslant n \leqslant N$ is as follows.
\begin{equation*}
\begin{aligned}
\frac{\partial ^\beta u(x,t)}{\partial x^{\beta}} \vert _{1} ^{n}&=\frac{(\Delta x)^{1-\beta}}{\Gamma(3-\beta)}[-a_{0} \frac{\partial u}{\partial x } |_{0}^{n}-\frac{a_{0}}{\Delta x} u_{1}^{n}+\frac{a_{0}}{\Delta x} u_{2}^{n}]+O(\Delta x)^{2-\beta}, \  (0\leqslant n \leqslant N),\\
\frac{\partial ^\beta u(x,t)}{\partial x^{\beta}} \vert _{j} ^{n}&=\frac{(\Delta x)^{1-\beta}}{\Gamma(3-\beta)}[-a_{j-1} \frac{\partial u}{\partial x } \vert_{0}^{n}+\frac{(a_{j-2}-a_{j-1})}{\Delta x} u_{1}^{n}\\
&+\Sigma_{k=2}^{j-1}\frac{(a_{j-k+1}-2a_{j-k}+a_{j-k-1})}{\Delta x} u_{k}^{n}+\frac{a_{1}-2a_{0}}{\Delta x} u_{j}^{n}+\frac{a_{0}}{\Delta x} u_{j+1}^{n}]\\
&+O(\Delta x)^{2-\beta}, \quad  (2 \leqslant j \leqslant J-1,\  0\leqslant n \leqslant N).
\end{aligned}
\end{equation*}
\end{lemma}
\begin{proof}
\begin{equation} \label{6}
\begin{aligned}
\frac{\partial ^\beta u(x,t)}{\partial x^{\beta}} \vert _{j} ^{n}=&\frac{1}{\Gamma (2- \beta)} \int_{0}^{x_{j}}{\frac{\partial ^2 u(s,t^{n})}{\partial s^{2}} (x_{j}-s)^{1-\beta}ds}\\
=&\frac{1}{\Gamma (2- \beta)} \Sigma_{k=1}^{j} [ (\frac{\frac{\partial u(x,t)}{\partial x} \vert_{k}^{n}-\frac{\partial u(x,t)}{\partial x} \vert_{k-1}^{n}}{\Delta x}\\
& \quad \quad \quad \quad \quad \quad \quad +O(\Delta x)) \int_{(k-1)\Delta x}^{k \Delta x}{(x_{j}-s)^{1-\beta}ds} ] \\
  =&\frac{(\Delta x)^{1-\beta}}{\Gamma (3-\beta)} \left\lbrace  -a_{j-1}\frac{\partial u}{\partial x}\vert_{0}^{n}+\Sigma_{k=1}^{j-1}{(a_{j-k}-a_{j-k-1})\frac{\partial u}{\partial x}\vert_{k}^{n}}+a_{0}\frac{\partial u}{\partial x}\vert_{j}^{n}\right\rbrace \\
 +&O(\Delta x)^{3-\beta}, \quad \quad \quad \quad  (1 \leqslant j \leqslant J-1, \quad 0 \leqslant n \leqslant N).
\end{aligned}
\end{equation}
Consider 
\begin{equation} \label{7}
\begin{aligned}
\frac{\partial u}{\partial x} \vert_{k}^{n}=\frac{u \vert _{k+1}^{n}-u \vert _{k}^{n}}{\Delta x}+O(\Delta x), \quad (1 \leqslant k \leqslant J-1, \quad 0 \leqslant n \leqslant N).
\end{aligned}
\end{equation}
Then, the relations \eqref{6} and \eqref{7} complete the proof.
\end{proof}\\

Assume
\begin{align} 
&\frac{\partial u}{\partial t} \vert_{j}^{n+\frac{1}{2}}=\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+O(\Delta t)^{2}, \ (0\leqslant n \leqslant N-1, \ 1 \leqslant j \leqslant J-1),  \label{8}\\
&\frac{\partial^{\beta} u}{\partial x^{\beta}} \vert_{j}^{n+\frac{1}{2}}=\frac{1}{2} \left[  \frac{\partial^{\beta} u}{\partial x^{\beta}} \vert_{j}^{n+1}+ \frac{\partial^{\beta} u}{\partial x^{\beta}} \vert_{j}^{n}   \right] +O(\Delta t)^{2}, \ (0\leqslant n \leqslant N-1, \ 1 \leqslant j \leqslant J-1). \label{9}
\end{align}
By disregarding the truncation errors, the discretization of equation \eqref{1} with conditions \eqref{2}-\eqref{4} at the grid point $ x_{j} \ (j=1, 2, ..., J-1)$ and time step $ (n+\frac{1}{2}) $ for $ 0 \leqslant n\leqslant N-1 $ using lemma \ref{lemma1} and relations \eqref{8} and \eqref{9} is as follows.
\begin{equation} \label{10}
\begin{aligned}
&u _{1}^{n+1}- \gamma_{1} \lbrace  \frac{ -a_{0}}{\Delta x} u_{1}^{n+1} + \frac{ a_{0}}{\Delta x} u_{2}^{n+1} \rbrace=u _{1}^{n}+\gamma_{1} \lbrace  \frac{ -a_{0}}{\Delta x} u_{1}^{n} + \frac{ a_{0}}{\Delta x} u_{2}^{n}\rbrace\\
&\quad  \quad +(\Delta t) f_{1}^{n+\frac{1}{2}}-\gamma_{1} a_{0}(\frac{\partial u}{\partial x} \vert_{0}^{n}+\frac{\partial u}{\partial x} \vert_{0}^{n+1}),\ 0 \leqslant n \leqslant N-1, \\
\end{aligned}
\end{equation}\\
\begin{equation} \label{11}
\begin{aligned}
&u _{j}^{n+1}- \gamma_{j} [  \frac{a_{j-2}-a_{j-1}}{\Delta x} u_{1}^{n+1}+\frac{\Sigma_{k=2}^{j-1}{(a_{j-k+1}-2a_{j-k}+a_{j-k-1})}}{\Delta x} u_{k}^{n+1}\\
& \quad \quad  \quad+\frac{a_{1}-2a_{0}}{\Delta x} u_{j}^{n+1}+\frac{a_{0}}{\Delta x} u_{j+1}^{n+1} ]  \\
&=u _{j}^{n}+ \gamma_{j} [ \frac{a_{j-2}-a_{j-1}}{\Delta x} u_{1}^{n}+\frac{\Sigma_{k=2}^{j-1}{(a_{j-k+1}-2a_{j-k}+a_{j-k-1})}}{\Delta x} u_{k}^{n}\\
& \quad \quad  \quad+\frac{a_{1}-2a_{0}}{\Delta x} u_{j}^{n}+\frac{a_{0}}{\Delta x} u_{j+1}^{n} ]+(\Delta t) f_{j}^{n+\frac{1}{2}}\\
&\quad \quad \quad -\gamma_{j} a_{j-1}(\frac{\partial u}{\partial x} \vert_{0}^{n}+\frac{\partial u}{\partial x} \vert_{0}^{n+1})  , \quad 2 \leqslant j\leqslant J-2,\ 0 \leqslant n \leqslant N-1,\\
\end{aligned}
\end{equation}\\
\begin{equation}\label{12}
\begin{aligned}
&{u^{n+1}_{{}_{J-1}}}- {\gamma_{{}_{J-1}}} [  \frac{{a_{{}_{{}_{J-3}}}}-{a_{{}_{J-2}}}}{\Delta x} u_{1}^{n+1}+\frac{\Sigma_{k=2}^{{}_{{}_{J-2}}}{[({a_{{}_{{}_{J-k}}}}-2{a_{{}_{J-k-1}}}+{a_{{}_{J-k-2}}})}]}{\Delta x} u_{k}^{n+1}\\
&\quad \quad  \quad +\frac{a_{1}-2a_{0}}{\Delta x} u_{{}_{J-1}}^{n+1} ]  \\
&={u^{n}_{{}_{J-1}}}+ {\gamma_{{}_{J-1}}} [ \frac{{a_{{}_{J-3}}}-{a_{{}_{J-2}}}}{\Delta x} u_{1}^{n}+\frac{\Sigma_{k=2}^{J-2}{[({a_{{}_{J-k}}}-2{a_{{}_{J-k-1}}}+{a_{{}_{J-k-2}})}}]}{\Delta x} u_{k}^{n}\\
&\quad  \quad +\frac{a_{1}-2a_{0}}{\Delta x} u_{{}_{J-1}}^{n}] +(\Delta t) f_{J-1}^{n+\frac{1}{2}}+{{\gamma}_{{}_{J-1}}} \frac{a_{0}}{\Delta x} (u_{{}_{J}}^{n}+u_{{}_{J}}^{n+1}) \\
&\quad \quad -{\gamma_{{}_{J-1}} }a_{{}_{{}_{J-2}}}(\frac{\partial u}{\partial x} \vert_{0}^{n}+\frac{\partial u}{\partial x} \vert_{0}^{n+1}), \quad 0 \leqslant n \leqslant N-1,\\
\end{aligned}
\end{equation}
where $\gamma_{j} = \frac{\alpha(x_{j}) (\Delta t) (\Delta x)^{1-\beta}}{2 \Gamma(3-\beta)}$ for $1 \leqslant  j \leqslant J-1 $, $ f_{j}^{n+\frac{1}{2}}=f(x_{j},t^{n+\frac{1}{2}}) $ for $ 1\leqslant j \leqslant J-1 $, and $ 0\leqslant n \leqslant N-1 $. Now the following theorem is easy to prove.\\
\begin{theorem} \label{consist}
The discretization of equation \eqref{1} with conditions \eqref{2}-\eqref{4} using lemma \ref{lemma1} and relations \eqref{8}-\eqref{9} is consistent with accuracy $ (O(\Delta x)^{2-\beta}+O(\Delta t)^{2}) $.
\end{theorem}

\section{Stability and convergence} \label{s-c}
The stability and convergence of schemes \eqref{10}-\eqref{12} are presented in this section. The idea to prove stability is taken from \cite{tadjeran2006}. Equations \eqref{10}-\eqref{12} can be considered as follows:
\begin{equation} \label{13}
(I-B)U^{n+1}=(I+B)U^{n}+F^{n+\frac{1}{2}}, \quad 0 \leqslant n \leqslant N-1,
\end{equation}
where 
\begin{equation*}
\begin{aligned}
&U^{n}=[u_{1}^{n}, u_{2}^{n}, ..., u_{J-1}^{n}]^{T}, \quad 0 \leqslant n \leqslant N-1,\\
&F^{n+\frac{1}{2}}=[(\Delta t){{f}_{1}^{n+\frac{1}{2}}}, (\Delta t){{f}_{2}^{n+\frac{1}{2}}}, ..., (\Delta t){f^{n+\frac{1}{2}}_{{}_{J-2}}}, (\Delta t){f^{n+\frac{1}{2}}_{{}_{J-1}}}+\frac{{{\gamma}_{J-1}} {{a}_{0}}}{\Delta x}({{u}^{n}_{{}_{J}}}+{{u}^{n+1}_{{}_{J}}})]^T\\
&\quad \quad \quad \quad -(\frac{\partial u }{\partial x} {{\vert}_{0}^{n}}+\frac{\partial u }{\partial x}{ {\vert}_{0}^{n+1}})[{\gamma_{1}}{ a_{0}}, {\gamma_{2}} {a_{1}}, ..., \gamma_{{}_{J-1}} a_{{}_{J-2}}]^T,\quad 0 \leqslant n \leqslant N-1,
\end{aligned}
\end{equation*}
$ I $ is a $ (J-1)\times (J-1) $ identity matrix, and in matrix $ B $, the elements $ B_{jk} \ (j, k=1, 2, ..., J-1) $ are as follows.
\begin{equation*}
\begin{aligned}
B_{jk}=
\begin{cases}
-\gamma_{1} \frac{a_{0}}{\Delta x}   & if \quad  k=j=1,\\
\gamma_{j} \frac{a_{1}-2a_{0}}{\Delta x}  &  if  \quad  k=j \ne 1,\\
\gamma_{j} \frac{a_{j-2}-a_{j-1}}{\Delta x}  & if \quad  k=1, 2 \leqslant j,\\
\gamma_{j} \frac{a_{j-k+1}-2a_{j-k}+a_{j-k-1}}{\Delta x} \quad & if \quad  2 \leqslant  k \leqslant  j-1, \ 2\leqslant j,\\
\gamma_{j} \frac{a_{0}}{\Delta x} & if \quad  k=j+1,\\
0 &  if \quad  j+1<k.
\end{cases}
\end{aligned}
\end{equation*}
\begin{theorem} \label{TS}
The finite difference discretization of equation \eqref{1} with conditions \eqref{2}-\eqref{4} defined by \eqref{10}-\eqref{12}  is unconditionally stable.
\end{theorem}
\begin{proof}
Equations  \eqref{10}-\eqref{12} are equivalent to \eqref{13}. First, it is argued that matrix $ B $ has the eigenvalues with a non--positive real--part. According to the Gershgorin theorem (\cite{datta2010} p. 294), for matrix $ B $, we have 
\begin{equation}\label{14}
\begin{cases}
\vert \lambda_{1}+\gamma_{1} \frac{a_{0}}{\Delta x} \vert \leqslant \gamma_{1} \frac{a_{0}}{\Delta x},\\
\vert \lambda_{j}+\gamma_{j} \frac{2a_{0}-a_{1}}{\Delta x} \vert \\
\quad \leqslant \gamma_{j} [\frac{a_{j-2}-a_{j-1}}{\Delta x}+\Sigma_{k=2}^{j-1} \vert {\frac{(a_{j-k+1}-2a_{j-k}+a_{j-k-1})}{\Delta x}} \vert+\frac{a_{0}}{\Delta x}], \ for  \ 2 \leqslant j \leqslant J-2,\\
\vert \lambda_{{}_{J-1}}+\gamma_{{}_{J-1}} \frac{2a_{0}-a_{1}}{\Delta x} \vert \leqslant \gamma_{{}_{J-1}} [\frac{a_{{}_{J-3}}-a_{{}_{J-2}}}{\Delta x}+\Sigma_{k=2}^{{}_{J-2}} \vert {\frac{(a_{{}_{J-k}}-2a_{{}_{J-k-1}}+a_{{}_{J-k-2}})}{\Delta x}} \vert],
\end{cases}
\end{equation}
where $ \lambda_{j} \ (1 \leqslant j \leqslant J-1) $ is the eigenvalue of the matrix $ B $. It is easy to show that $ a_{0} > a_{1} > ... > a_{n} $ and $ lim_{n\rightarrow \infty} a_{n}=0$. Also, we can show that
\begin{equation*}
 (a_{n}-a_{n+1} )> (a_{n+1}-a_{n+2} ),  \quad \quad   n=1, 2, ... .
\end{equation*}
Therefore, for $ 2 \leqslant j \leqslant J-1 $ and $ k=2, 3, ..., j-1 $, we have 
\begin{equation*}
(a_{{}_{j-k+1}}-2a_{{}_{j-k}}+a_{{}_{j-k-1}}) >0 . 
\end{equation*}
Therefore,
\begin{equation*}
\begin{aligned}
 \Sigma_{k=2}^{{}_{j-1}}  {\vert(a_{{}_{j-k+1}}-2a_{{}_{j-k}}+a_{{}_{j-k-1}})\vert} &= \Sigma_{k=2}^{{}_{j-1}} ((a_{{}_{j-k+1}}-a_{{}_{j-k}})-(a_{{}_{j-k}}-a_{{}_{j-k-1}}))\\
 &=(a_{j-1}-a_{j-2})+(a_{0}-a_{1}).
\end{aligned}
\end{equation*}
 Now the relations \eqref{14} can be written as follows.
\begin{equation}\label{15}
\begin{aligned}
\begin{cases}
\vert \lambda_{1}+\gamma_{1} \frac{a_{0}}{\Delta x} \vert \leqslant \gamma_{1} \frac{a_{0}}{\Delta x},\\
\vert \lambda_{j}+\gamma_{j} \frac{2a_{0}-a_{1}}{\Delta x} \vert \leqslant \gamma_{j} \frac{2a_{0}-a_{1}}{\Delta x}, \quad 2\leqslant j \leqslant J-2,\\
\vert \lambda_{J-1}+\gamma_{{}_{J-1}} \frac{2a_{0}-a_{1}}{\Delta x} \vert \leqslant \gamma_{{}_{J-1}} \frac{a_{0}-a_{1}}{\Delta x}.
\end{cases}
\end{aligned}
\end{equation}
It is obvious that matrix $ B $ is invertible, so matrix $ B $ has the non--zero eigenvalues. Therefore according to \eqref{15}, matrix $ B $ has the eigenvalues with the non--positive real--part.\\
Now, $ \lambda_{j} $ is an eigenvalue of the matrix $ B $ if and only if $ \frac{1+\lambda_{j}}{1-\lambda_{j}} $ is an eigenvalue of the matrix $ (I-B)^{-1}(I+B) $. Since the real part of $ \lambda_{j} $ is not positive, $ \vert \frac{1+\lambda_{j}}{1-\lambda_{j}} \vert <1$. Thus, the system of equation \eqref{13} is unconditionally stable.
\end{proof}\\

By using Lax's equivalence theorem \cite{laxandrichtmyer1956}, theorem \ref{consist} and theorem \ref{TS} indicate that our proposed scheme \eqref{13} is convergent.
\section{Numerical tests} \label{test}
Some numerical tests are presented in this section, to check the validity of our proposed scheme.  We measure the accuracy of the proposed method by assuming $\Delta x$ and $\Delta t$, using the following maximum absolute error 
\begin{equation*}
L_{\infty}(\Delta x, \Delta t)=max_{{}_{{1 \leqslant j \leqslant J-1},{1 \leqslant n \leqslant N} }}{\vert \widehat{u}_{j}^{n}-u_{j}^{n} \vert},
\end{equation*}
where $\widehat{u}_{j}^{n}$, and $  u_{j}^{n}$ are the approximation and the exact solutions of equation \eqref{1} with conditions \eqref{2}--\eqref{4} at $ x_{j} $ and time $ t ^{n}$, respectively. To test our proposed method, we consider the following three examples in which the exact solutions are available.
\begin{example} \label{example1}
Assume the equation \cite{khader2011}
\begin{equation*}
\frac{\partial u(x,t)}{\partial t}=\alpha(x) \frac{\partial^{1.8} u(x,t)}{\partial x^{1.8}}+f(x,t), \quad 0<x<1, \quad 0<t\leqslant 1,
\end{equation*}
where $ \alpha(x)=\Gamma(1.2) x^{1.8}$, $ f(x,t)=(6x^3-3x^2)e^{-t} $. The exact solution $ u(x,t)=(x^{2}-x^{3})e^{-t} $ is used to consider conditions \eqref{2}--\eqref{4}.
\end{example}
\begin{example} \label{example2}
Assume the equation \cite{saadatmandi2011}
\begin{equation*}
\frac{\partial u(x,t)}{\partial t}=\alpha(x) \frac{\partial^{1.5} u(x,t)}{\partial x^{1.5}}+f(x,t), \quad 0<x<1, \quad 0<t\leqslant 1,
\end{equation*}
where $ \alpha(x)=\Gamma(1.5) x^{0.5}$, $ f(x,t)=(x^2+1) cos(t+1)-2x sin(t+1)$. The exact solution $ u(x,t)=(x^2+1) sin(t+1) $  is used to consider conditions \eqref{2}--\eqref{4}.
\end{example}
\begin{example} \label{example3}
Assume the equation \cite{tadjeran2006}
\begin{equation*}
\frac{\partial u(x,t)}{\partial t}=\alpha(x) \frac{\partial^{1.8} u(x,t)}{\partial x^{1.8}}+f(x,t), \quad 0<x<1, \quad 0<t\leqslant 1,
\end{equation*}
where $ \alpha(x)=\Gamma(2.2) \frac{x^{2.8}}{6}$, $ f(x,t)=-(1+x)e^{-t}x^3$. The exact solution $ u(x,t)=e^{-t}x^3$  is used to consider conditions \eqref{2}--\eqref{4}.
\end{example}
The results of our proposed method for Examples \ref{example1}, \ref{example2}, and \ref{example3} are presented in Tables \ref{Example1}, \ref{Example2}, and \ref{Example3}, respectively. The second columns of these Tables show that the maximum absolute error,  with  different values $\Delta t $ and $ \Delta x $, is small enough and reduces as the grids are refined.

To test the rate of convergence of our proposed scheme, we started with $ \Delta x=\Delta t=\frac{1}{10} $ and obtained numerical solutions for Examples \ref{example1}, \ref{example2}, and \ref{example3}. We then repeated the computations using finer grids. Here, the Error rate is defined by the ratio of the errors as refining the grids, as follows.
\begin{equation*}
Error \ rate=\frac{L_{\infty}((\Delta x)_{1}, (\Delta t)_{1})}{L_{\infty}((\Delta x)_{_{2}}, (\Delta t)_{_{2}})},
\end{equation*}
where$ (\Delta x)_{_{2}}= (\Delta t)_{_{2}}<(\Delta x)_{_{1}}=(\Delta t)_{_{1}}$. According to the third columns of Tables \ref{Example1}, \ref{Example2}, and \ref{Example3}, the behavior of errors is (almost) linear. It means that the ratio of $ L_{\infty}((\Delta x)_{_{1}}, (\Delta t)_{_{1}})$ to $L_{\infty}((\Delta x)_{_{2}}, (\Delta t)_{_{2}})$ is approximately equal to the ratio of $(\Delta x)_{_{1}}$ or $(\Delta t)_{_{_{1}}}$ to $(\Delta x)_{_{2}}$ or $(\Delta t)_{_{2}}$, where $(\Delta x)_{_{2}}= (\Delta t)_{_{2}}<(\Delta x)_{_{1}}= (\Delta t)_{_{1}}$. 
\begin{table} [h!]
\centering
% table caption is above the table
\caption{The maximum absolute errors and error rates  with  different values $\Delta t $ and $ \Delta x $  in Example 1. }\label{Example1}
% For LaTeX tables use
\begin{tabular}{lll}
\hline\noalign{\smallskip}
$\Delta t ,\Delta x$ &$L_{\infty}(\Delta x, \Delta t) $  & Error rate \\
\noalign{\smallskip}\hline\noalign{\smallskip}
 $\frac{1}{10}$ & $ 7.7904 e-3 $ & $ -$\\  
$  $ &$  $ &$  $ \\
 $\frac{1}{20}$ & $4.6223 e-3  $ & $ 1.68 \thickapprox \frac{20}{10}$\\  
$  $ &$  $ &$  $ \\
 $\frac{1}{50}$ & $ 2.1667 e-3 $ & $2.13 \thickapprox \frac{50}{20}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{100}$ & $ 1.1842 e-3 $ & $1.83 \thickapprox \frac{100}{50} $\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{200}$ & $6.3542 e-4  $ & $ 1.86 \thickapprox \frac{200}{100}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{500}$ & $ 2.7344 e-4 $ & $2.32 \thickapprox \frac{500}{200}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{1000}$ & $ 1.4287 e-4 $ & $1.91 \thickapprox \frac{1000}{500} $\\ 
$  $ &$  $ &$  $ \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table} 
\begin{table}[h!]
\centering
% table caption is above the table
\caption{The maximum absolute errors and error rates  with  different values $\Delta t $ and $ \Delta x $  in Example 2. }\label{Example2}
% For LaTeX tables use
\begin{tabular}{lll}
\hline\noalign{\smallskip}
$\Delta t ,\Delta x$ &$L_{\infty}(\Delta x, \Delta t) $  & Error rate \\
\noalign{\smallskip}\hline\noalign{\smallskip}
 $\frac{1}{10}$ & $ 5.4187 e-2 $ & $ -$\\  
$  $ &$  $ &$  $ \\
 $\frac{1}{20}$ & $  2.6464 e-2$ & $ 2.05 \thickapprox \frac{20}{10}$\\  
$  $ &$  $ &$  $ \\
 $\frac{1}{50}$ & $  1.0139 e-2$ & $ 2.61\thickapprox \frac{50}{20}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{100}$ & $ 4.9237 e-3 $ & $ 2.06\thickapprox \frac{100}{50} $\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{200}$ & $2.4043 e-3 $ & $2.05 \thickapprox \frac{200}{100}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{500}$ & $ 9.3982 e-4$ & $ 2.56\thickapprox \frac{500}{200}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{1000}$ & $4.6414 e-4 $ & $2.04\thickapprox \frac{1000}{500} $\\ 
$  $ &$  $ &$  $ \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table} 
\begin{table}[h!]
\centering
% table caption is above the table
\caption{The maximum absolute errors and error rates  with  different values $\Delta t $ and $ \Delta x $  in Example 3. }\label{Example3}
% For LaTeX tables use
\begin{tabular}{lll}
\hline\noalign{\smallskip}
$\Delta t ,\Delta x$ &$L_{\infty}(\Delta x, \Delta t) $  & Error rate \\
\noalign{\smallskip}\hline\noalign{\smallskip}
$\frac{1}{10}$ & $  4.3914 e-3$ & $ -$\\  
$  $ &$  $ &$  $ \\
 $\frac{1}{20}$ & $ 2.4129 e-3 $ & $  1.82\thickapprox \frac{20}{10}$\\  
$  $ &$  $ &$  $ \\
 $\frac{1}{50}$ & $  1.1076 e-3$ & $2.18 \thickapprox \frac{50}{20}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{100}$ & $ 5.7498 e-4 $ & $ 1.93 \thickapprox \frac{100}{50} $\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{200}$ & $ 3.0367 e-4$ & $ 1.89 \thickapprox \frac{200}{100}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{500}$ & $1.2873 e-4 $ & $ 2.36 \thickapprox \frac{500}{200}$\\ 
$  $ &$  $ &$  $ \\
 $\frac{1}{1000}$ & $6.6693 e-5 $ & $ 1.93\thickapprox \frac{1000}{500} $\\ 
$  $ &$  $ &$  $ \\
\noalign{\smallskip}\hline
\end{tabular}
\end{table} 

\section{Conclusion} \label{conclusion}
In this paper, a numerical solution for a one-dimensional space fractional diffusion equation in the sense of Caputo with the initial condition, Neumann, and Dirichlet boundary conditions was presented using the finite difference method. The consistency and stability of the proposed method were proven. Then, the convergence of this scheme was obtained using Lax's equivalence theorem. Numerical tests show that the proposed scheme is accurate enough, and the behavior of errors is nearly linear.

\newpage
 \begin{center}
 \begin{thebibliography}{99} 
\bibitem{agarwal2018} P. Agarwal and A. A. El-Sayed, Non-standard finite difference and Chebyshev collocation methods for solving fractional diffusion equation, {\it Physica A: Statistical Mechanics and its Applications}, 500 (2018), 40-49. 
\bibitem{datta2010} B. N. Datta, {\it Numerical Linear Algebra and Applications}, 2nd edition, SIAM, (2010). \bibitem{khader2011} M. M. Khader, On the numerical solutions for the fractional diffusion equation, {\it Communications in Nonlinear Science and Numerical Simulation}, 16(6): (2011), 2535-2542. 
\bibitem{kheybari2019} S. Kheybari, M. T. Darvishi, and M. S. Hashemi, Numerical simulation for the space-fractional diffusion equations, {\it Applied Mathematics and Computation}, 348 (2019), 57-69. 
\bibitem{laxandrichtmyer1956} P. D. Lax and R. D. Richtmyer, Survey of the stability of linear finite difference equations, {\it Communications on Pure and Applied Mathematics}, 9(2): (1956), 267-293.
\bibitem{meerschaert2006} M. M. Meerschaert, H-P. Scheffler, and C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation, {\it Journal of Computational Physics}, 211(1): (2006), 249-261. 
\bibitem{ren2013} R-f. Ren, H-b. Li, W. Jiang, and M-y. Song, An efficient Chebyshev-tau method for solving the space fractional diffusion equations, {\it Applied Mathematics and Computation}, 224: (2013), 259-267. 
 \bibitem{saadatmandi2011} A. Saadatmandi and M. Dehghan, A tau approach for solution of the space fractional diffusion equation, {\it Computers and Mathematics with Applications}, 62(3): (2011), 1135-1142.
 \bibitem{safdari2019} H. Safdari, H. Mesgarani, M. Javidi, and Y. E. Aghdam, Convergence analysis of the space fractional-order diffusion equation based on the compact finite difference scheme, {\it Computational and Applied Mathematics}, 39(2): (2020), 1-15.
\bibitem{tadjeran2006} C. Tadjeran, M. M. Meerschaert, and H-P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, {\it Journal of Computational Physics}, 213(1): (2006), 205-213. \end{thebibliography} \end{center}



{\small

\noindent{\bf First Author (Haniye Hajinezhad) }

\noindent Department of Mathematics

\noindent Assistant Professor of Mathematics

\noindent Payame Noor University,

\noindent Tehran, Iran

\noindent E-mail: H.Hajinezhad@pnu.ac.ir}\\

{\small
\noindent{\bf  Second Author (Ali Reza Soheili)  }

\noindent  Department of Mathematics

\noindent Professor of Mathematics

\noindent  Faculty of Mathematical Sciences, Ferdowsi University of Mashhad

\noindent  Mashhad, Iran

\noindent E-mail: soheili@um.ac.ir}\\



\end{document}