## 4 An Implicit Finite Difference Method
As usual, we supposed that the solution was sufficiently smooth and discretized its first \(\partial/\partial x\) and second order \(\partial^{2}/\partial x^{2}\) spatial derivatives by means of the second order, three-point central difference formula so that \[\frac{\partial}{\partial x}u(x_{j},t^{n})=\frac{u(x_{j+1},t^{n})-u(x_{j-1},t^{n})}{ 2\Delta x}+\mathcal{O}(\Delta x^{2}) \tag{6}\]  
\[\frac{\partial^{2}}{\partial x^{2}}u(x_{j},t^{n})=\frac{u(x_{j+1},t^{n})-2u(x_{ j},t^{n})+u(x_{j-1},t^{n})}{\Delta x^{2}}+\mathcal{O}(\Delta x^{2}). \tag{7}\]  
According to the discretization for the Caputo derivative (4), we approximated the time fractional derivative in (1) as follows  
\[\frac{\partial^{\alpha}}{\partial t^{\alpha}}u(x_{j},t^{n})= \tag{8}\] \[\frac{1}{\Gamma(2-\alpha)}\sum_{k=1}^{n}\frac{u(x_{j},t^{k})-u(x_{ j},t^{k-1})}{t^{k}-t^{k-1}}\left[(t^{n}-t^{k-1})^{1-\alpha}-(t^{n}-t^{k})^{1- \alpha}\right]+R_{j}^{n},\]  
where  
\[\frac{1}{\Gamma(1-\alpha)}\int_{t^{k-1}}^{t^{k}}(t^{n}-s)^{-\alpha}ds=\frac{ 1}{\Gamma(2-\alpha)}[(t^{n}-t^{k-1})^{1-\alpha}-(t^{n}-t^{k})^{1-\alpha}].\]  
Replacing \(u(x_{j},t^{n})\) with its numerical approximation, \(U_{j}^{n}\), and neglecting the local truncation errors, the time-fractional advection-diffusion Equation (1) was discretized as follows  
\[\frac{1}{\Gamma(2-\alpha)}\sum_{k=1}^{n}T_{n,k}(U_{j}^{k}-U_{j}^{k-1})+K_{1} \frac{U_{j+1}^{n}-U_{j-1}^{n}}{2\Delta x}-K_{2}\frac{U_{j+1}^{n}-2U_{j}^{n}+U_{ j-1}^{n}}{\Delta x^{2}}=f_{j}^{n}, \tag{9}\]  
for \(1\leq n\leq N\) and \(1\leq j\leq J-1\), where \(f_{j}^{n}=f(x_{j},t^{n})\) is the source term and  
\[T_{n,k}=\frac{(t^{n}-t^{k-1})^{1-\alpha}-(t^{n}-t^{k})^{1-\alpha}}{t^{k}-t^{k -1}}\qquad 1\leq k\leq n,\qquad 1\leq n\leq N.\]  
The initial and boundary conditions can be rewritten as  
\[U_{j}^{0}=\phi(x_{j})\qquad 0\leq j\leq J\] \[U_{0}^{n}=\phi(t^{n})\qquad U_{J}^{n}=\psi(t^{n})\qquad 0\leq n\leq N. \tag{10}\]  
For any temporal meshes on \([0,T]\) and for any \(n=1,2,\cdots,N\), it holds that [22]  
\[T_{n,k}>0\qquad\qquad T_{n,k}>T_{n,k-1},\qquad 1\leq k\leq n.\]  
Taking into account that \(T_{n,n}=(\Delta t^{n})^{-\alpha}\), we can write  
\[\frac{\partial^{2}}{\partial x^{2}}u(x_{j},t^{n})=\frac{u(x_{j+1},t^{n})-2u(x_{ j},t^{n})+u(x_{j-1},t^{n})}{\Delta x^{2}}+\mathcal{O}(\Delta x^{2}). \tag{7}\]0  
Then, the following implicit finite difference scheme was obtained:  
\[\frac{\partial^{2}}{\partial x^{2}}u(x_{j},t^{n})=\frac{u(x_{j+1},t^{n})-2u(x_{ j},t^{n})+u(x_{j-1},t^{n})}{\Delta x^{2}}+\mathcal{O}(\Delta x^{2}). \tag{7}\]1 \[\frac{\partial^{2}}{\partial x^{2}}u(x_{j},t^{n})=\frac{u(x_{j+1},t^{n})-2u(x_{ j},t^{n})+u(x_{j-1},t^{n})}{\Delta x^{2}}+\mathcal{O}(\Delta x^{2}). \tag{7}\]2  
where we set  
\[\frac{\partial^{2}}{\partial x^{2}}u(x_{j},t^{n})=\frac{u(x_{j+1},t^{n})-2u(x_{ j},t^{n})+u(x_{j-1},t^{n})}{\Delta x^{2}}+\mathcal{O}(\Delta x^{2}). \tag{7}\]3The Equation (11) can be written in vectorial form as  
\[\frac{\partial^{2}}{\partial x^{2}}u(x_{j},t^{n})=\frac{u(x_{j+1},t^{n})-2u(x_{ j},t^{n})+u(x_{j-1},t^{n})}{\Delta x^{2}}+\mathcal{O}(\Delta x^{2}). \tag{7}\]4  
where \(K\) is a tridiagonal matrix and where \(\mathcal{L}\) denotes the following difference operator  
\[\frac{\partial^{2}}{\partial x^{2}}u(x_{j},t^{n})=\frac{u(x_{j+1},t^{n})-2u(x_{ j},t^{n})+u(x_{j-1},t^{n})}{\Delta x^{2}}+\mathcal{O}(\Delta x^{2}). \tag{7}\]5  
Here, and in the following text, we assume the convention that the summation is equal to zero if the lower bound is larger that the upper bound. The obtained method (12) is implicit. In order to compute the numerical solution, \(U^{n}\), a system with the tridiagonal coefficients matrix, \(K\), has to be solved. It is interesting to note that the operator, \(\mathcal{L}\), is a kind operator with memory, due to the non-local character of the fractional derivative. This means that the effect on \(U\) at time \(t^{n}\), \(U^{n}\) depends on all the previous values, \(U^{0},U^{1},\cdots,U^{n-1}\), evaluated at all the previous time, \(t^{0},t^{1},\cdots,t^{n-1}\).