Courtesy of Santos Bravo Yuste  ☺   
Dpt. Fisica. UEx. Badajoz. Spain  
email:  santos@unex.es
http://www1.unex.es/eweb/fisteor/santos/sby.html

Weighted average finite difference methods for fractional diffusion equations

The fractional version of the  Crank-Nicolson method is just a particular case of the weighted average method in which the weight parameter λ  is  1/2

Here I provide a   Mathematica notebook that uses the weighted-averaged difference method  for solving some class of fractional diffusion equations. The method is discussed in 
	Ref[1] = S. B. Yuste,  Weighted average finite difference methods for fractional diffusion equations, Proceedings of FDA'04, pages 335-340 (2004)  [cs.NA/0408053] 
and
	Ref [2] = S. B. Yuste, Weighted average finite difference methods for fractional diffusion equations, Journal of Computational Physics 216 (2006) 264-274

The fractional diffusion equation (FDE) we are going to solve is  (see Eq. (1) of Ref[2]) : 

	(1)		((∂u)/(∂t))= K  (((Subscript[, o])D)^(1-γ)) (((∂^2)u)/(∂(x^2)))	     
	
In particular, as an example, we are going to use the initial condition u(x,0)=x(1-x)  and the  Dirichlet boundary conditions u(0,t)=u(1,t)=0

The solution u evaluated at position (Subscript[x, j])  at time (Subscript[t, m]),   u((Subscript[x, j]),(Subscript[t, m])), is denoted here as u[j,m,γ], where the third index γ  is the anomalous diffusion exponent of the FDE.

The exact solution

uExactx1mx[x,t,γ,Kg,numterm] is the  exact solution  of the FDE (1) written as a (finite) series with numterm terms (the exact solution corresponds to numterm→∞)

uExactx1mx[x_,t_,γ_,Kg_,numterm_]:=8/Pi^3*Sum[1/(2*n+1)^3*Sin[(2*n+1)*Pi*x]*MittagLefflerE[γ,-Kg*(2*n+1)^2*Pi^2*t^γ ],{n,0,numterm}]

(The  Mittag-Leffler function for γ = 1/2  can be written in a simpler way:   MittagLefflerE[1/2,z]=(e^(z^2)) Erfc[-z])

The Weighted averaged (WA) numerical algorithm  (see Eq.(19) of Ref[2]):

BDF1 Coefficients:  w[k_, alfa_]≡(Subscript[ω, k]^((α)))  (see formula  (13) of Ref[2]) :

w[0,alfa_]:=1;w[k_,alfa_]:=w[k,alfa]=(1-(alfa+1)/k)*w[k-1,alfa]

In what follows, 
	lambda≡λ  (weight factor; see Eq. (17) of [2]) ,  
	 sfrac≡S  (parameter defined in Eq. (19b) of  [2]) and 
	 Nx determines the space discretization of the interval in which the equation is going to be solve (the total number of points (Subscript[x, j]) is   2Nx+1)

RR[j_,m_,gamma_]:=Module[{kk},
kk=Sum[((1-lambda)*w[k+1,1-gamma]+lambda*w[k,1-gamma])*(uWA[j-1,m-1-k,gamma]-2*uWA[j,m-1-k,gamma]+uWA[j+1,m-1-k,gamma]),{k,0,m-1}];
uWA[j,m-1,gamma]+sfrac*kk]

uWA[j, m, gamma]  provides the value of the numeral solution at the point (Subscript[x, j])  at time  (Subscript[t, m])  when the  anomalous diffusion exponent γ is gamma.

uWALista[m_,gamma_]:=uWALista[m,gamma]=LinearSolve[SparseArray[{Band[{1,2}]->Table[Which[j==Nx-1,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}],Band[{1,1}]->Table[Which[j==-Nx,1,j==Nx,1,True,(1+2*sfrac*w[0,1-gamma]*(1-lambda))],{j,-Nx,Nx}],Band[{2,1}]->Table[Which[j==-Nx,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}]},2*Nx+1],Table[Which[j==-Nx,uL[m*ht],j==Nx,uR[m*ht],True,RR[j,m,gamma]],{j,-Nx,Nx}]]

uWA[j_,m_,gamma_]:=uWA[j,m,gamma]=uWALista[m,gamma][[Nx+1+j]]

soluWA[t, gamma]   is a function that generates a table  with the solution  u  evaluated at time t. The table has the format {x, u (x, t)}

soluWA[t_,gamma_]:=Module[{mfin,kk,uMEfinalTab,eps=$MachineEpsilon},
mfin=Floor[(t+eps)/ht];kk=Timing[Do[uWAfinalTab=Table[{(Nx+j)*deltax,uWA[j,m,gamma]},{j,-Nx,Nx}],{m,0,mfin}]][[1]];Print["t=",mfin*ht," : mfin=",mfin," : ht=",ht," : Delta x=",deltax," : Weight factor λ=",lambda," : Computing time=",kk," : Weighted-Averaged method"];uWAfinalTab]

Examples

Numerical solution of the FDE by means of the weighted-averaged method with weight factor λ=1/2, i.e.,  by means of the fractional Crank-Nicolson method,  for 
	γ=1/2,   K=1,   ,   sfrac=0.7  ,   Δx=1/20 , Δt=6.25×(10^(-6)).
The fractional Crank-Nicolson method is always stable

Here we clear the saved values of uWALista and uWA  (the evaluation of this cell is required if you are  going to work with several examples in a Mathematica session)

Clear[uWALista, uWA];
uWALista[m_,gamma_]:=uWALista[m,gamma]=LinearSolve[SparseArray[{Band[{1,2}]->Table[Which[j==Nx-1,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}],Band[{1,1}]->Table[Which[j==-Nx,1,j==Nx,1,True,(1+2*sfrac*w[0,1-gamma]*(1-lambda))],{j,-Nx,Nx}],Band[{2,1}]->Table[Which[j==-Nx,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}]},2*Nx+1],Table[Which[j==-Nx,uL[m*ht],j==Nx,uR[m*ht],True,RR[j,m,gamma]],{j,-Nx,Nx}]]
uWA[j_,m_,gamma_]:=uWA[j,m,gamma]=uWALista[m,gamma][[Nx+1+j]]

lambda=weight factor=λ

lambda=0.5;

2*Nx+1=number of spatial points

Nx=10;

Parameters of this example:  gamma = 1/2, Kg = 1, sfrac = 1, deltax = 1./(2*Nx)   
(Notation:   gamma=γ,  ht=Δt  ,   Kg=K,   sfrac=S,  deltax=Δx )

{gamma=0.5,deltax=1./(2*Nx),Kg=1.;sfrac=1,ht=(sfrac*deltax^2/Kg)^(1/gamma)}

According to Eq. (27) of Ref[2], the critical stability value for sfrac≡S=K(((Δt))^γ)/(((Δx))^2) is given by 
                          (Subscript[S, × ])= 1/(2*(2*lambda-1)*  w(-1,1-gamma)
As we use the BDF1 coefficients, we have
                                    w(-1,1-gamma)=2^(1-gamma)

The present WA numerical method is stable if   
		1/sfracEfect ≥  1/(Subscript[S, × ])
with  
		sfracEfect=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2   
(see Eq. (25) of Ref[2])

{(Sin[(2*Nx-1)*Pi/(4*Nx)])^2//N;sfracEfectivo=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2;scritico=1/(2*(2*lambda-1)*2^(1-gamma));1/sfracEfectivo,1/scritico}

Here we see that   1/sfracEfect=((Sec[(π/40)])^2) is larger than  1/(Subscript[S, ×])=0.  Therefore, the algorithm is stable

We set the initial condition : u(x,t=0)= x(1-x)

uWA[nx_,0,alfa_]:=(nx+Nx)*deltax*(1-(nx+Nx)*deltax)

Table[{nx,uWA[nx,0,1]},{nx,-Nx,Nx}]//N

Boundary condition: u(0,t)=u(1,t)=0.

uL[t_]:=0;
uR[t_]:=0;

uWA[-Nx,m_,alfa_]:=uL[m*ht];
uWA[Nx,m_,alfa_]:=uR[m*ht];

The solution at the points x_j after mfin= 2 timesteps.

mfin=2;uWALista[mfin,gamma]

The solution in the format   {....{(Subscript[x, j]), u((Subscript[x, j]), (Subscript[t, m]))}, {(Subscript[x, j+1]), u((Subscript[x, j+1]), (Subscript[t, m]))},....}   where m≡mfin

mfin=1;soluWA[mfin*ht,gamma]

mfin=10;soluWA[mfin*ht,gamma]

mfin=100;soluWA[mfin*ht,gamma]

mfin=200;figNum=ListPlot[soluWA[mfin*ht,gamma],PlotRange->All,PlotStyle->PointSize[0.02]]

tv=200*ht;figExact=Plot[uExactx1mx[x,tv,1/2,1,20],{x,0,1},PlotStyle->{Red,Thick},PlotRange->All]

Show[figNum,figExact,PlotRange->All]

Numerical solution of the FDE by means of the weighted-averaged method with weight factor λ=1/2, i.e.,  by means of the fractional Crank-Nicolson method,  for 
	γ=3/4,   K=1,   ,   sfrac=0.7  ,   Δx=1/20 , Δt=6.25×(10^(-6)).
The fractional Crank-Nicolson method is always stable

Here we clear the saved values of uWALista and uWA  (the evaluation of this cell is required if you are  going to work with several examples in a Mathematica session)

Clear[uWALista, uWA];
uWALista[m_,gamma_]:=uWALista[m,gamma]=LinearSolve[SparseArray[{Band[{1,2}]->Table[Which[j==Nx-1,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}],Band[{1,1}]->Table[Which[j==-Nx,1,j==Nx,1,True,(1+2*sfrac*w[0,1-gamma]*(1-lambda))],{j,-Nx,Nx}],Band[{2,1}]->Table[Which[j==-Nx,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}]},2*Nx+1],Table[Which[j==-Nx,uL[m*ht],j==Nx,uR[m*ht],True,RR[j,m,gamma]],{j,-Nx,Nx}]]
uWA[j_,m_,gamma_]:=uWA[j,m,gamma]=uWALista[m,gamma][[Nx+1+j]]

lambda=weight factor=λ

lambda=0.5;

2*Nx+1=number of spatial points

Nx=10;

Parameters of this example:  gamma = 3/4, Kg = 1, sfrac = 1, deltax = 1./(2*Nx)   
(Notation:   gamma=γ,  ht=Δt  ,   Kg=K,   sfrac=S,  deltax=Δx )

{gamma=0.75,deltax=1./(2*Nx),Kg=1.;sfrac=1,ht=(sfrac*deltax^2/Kg)^(1/gamma)}

According to Eq. (27) of Ref[2], the critical stability value for sfrac≡S=K(((Δt))^γ)/(((Δx))^2) is given by 
                          (Subscript[S, × ])= 1/(2*(2*lambda-1)*  w(-1,1-gamma)
As we use the BDF1 coefficients, we have
                                    w(-1,1-gamma)=2^(1-gamma)

The present WA numerical method is stable if   
		1/sfracEfect ≥  1/(Subscript[S, × ])
with  
		sfracEfect=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2   
(see Eq. (25) of Ref[2])

{(Sin[(2*Nx-1)*Pi/(4*Nx)])^2//N;sfracEfectivo=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2;scritico=1/(2*(2*lambda-1)*2^(1-gamma));1/sfracEfectivo,1/scritico}

Here we see that   1/sfracEfect=((Sec[(π/40)])^2) is larger than  1/(Subscript[S, ×])=0.  Therefore, the algorithm is stable

We set the initial condition : u(x,t=0)= x(1-x)

uWA[nx_,0,alfa_]:=(nx+Nx)*deltax*(1-(nx+Nx)*deltax)

Table[{nx,uWA[nx,0,1]},{nx,-Nx,Nx}]//N

Boundary condition: u(0,t)=u(1,t)=0.

uL[t_]:=0;
uR[t_]:=0;

uWA[-Nx,m_,alfa_]:=uL[m*ht];
uWA[Nx,m_,alfa_]:=uR[m*ht];

The solution at the points x_j after mfin= 2 timesteps.

mfin=2;uWALista[mfin,gamma]

The solution in the format   {....{(Subscript[x, j]), u((Subscript[x, j]), (Subscript[t, m]))}, {(Subscript[x, j+1]), u((Subscript[x, j+1]), (Subscript[t, m]))},....}   where m≡mfin

mfin=1;soluWA[mfin*ht,gamma]

mfin=10;soluWA[mfin*ht,gamma]

mfin=100;soluWA[mfin*ht,gamma]

mfin=200;figNum=ListPlot[soluWA[mfin*ht,gamma],PlotRange->All,PlotStyle->PointSize[0.02]]

tv=200*ht;figExact=Plot[uExactx1mx[x,tv,3/4,1,20],{x,0,1},PlotStyle->{Red,Thick},PlotRange->All]

Show[figNum,figExact,PlotRange->All]

Numerical solution of the FDE by means of the weighted-averaged method with weight factor λ=0.8  for 
	γ=1/2,   K=1,   ,   sfrac=0.7  ,   Δx=1/20 , Δt=3.0625000000000007`*^-6  
The algorithm is unstable in this case

Here we clear the saved values of uWALista and uWA  (the evaluation of this cell is required if you are  going to work with several examples in a Mathematica session)

Clear[uWALista, uWA];
uWALista[m_,gamma_]:=uWALista[m,gamma]=LinearSolve[SparseArray[{Band[{1,2}]->Table[Which[j==Nx-1,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}],Band[{1,1}]->Table[Which[j==-Nx,1,j==Nx,1,True,(1+2*sfrac*w[0,1-gamma]*(1-lambda))],{j,-Nx,Nx}],Band[{2,1}]->Table[Which[j==-Nx,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}]},2*Nx+1],Table[Which[j==-Nx,uL[m*ht],j==Nx,uR[m*ht],True,RR[j,m,gamma]],{j,-Nx,Nx}]]
uWA[j_,m_,gamma_]:=uWA[j,m,gamma]=uWALista[m,gamma][[Nx+1+j]]

lambda=weight factor=λ

lambda=0.8;

2*Nx+1=number of spatial points

Nx=10;

Parameters of the example :  gamma = 1/2, Kg = 1, sfrac = 0.7, deltax = 1./(2*Nx)   
Notation:   gamma=γ,  ht=Δt  ,   Kg=K  ,   sfrac=(Subscript[S,  ]) ,  deltax=Δx

{gamma=0.5,deltax=1./(2*Nx),Kg=1.;sfrac=0.7,ht=(sfrac*deltax^2/Kg)^(1/gamma)}

According to Eq. (27) of Ref[2], the critical stability value for sfrac≡S=K(((Δt))^γ)/(((Δx))^2) is given by 
                          (Subscript[S, × ])= 1/(2*(2*lambda-1)*  w(-1,1-gamma)
As we use the BDF1 coefficients, we have
                                    w(-1,1-gamma)=2^(1-gamma)

The present numerical method is stable if   1/sfracEfect ≥  1/(Subscript[S, × ])with  sfracEfect=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2   
(see Eq. (25) of Ref[2])

{(Sin[(2*Nx-1)*Pi/(4*Nx)])^2//N;sfracEfectivo=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2;scritico=1/(2*(2*lambda-1)*2^(1-gamma));1/sfracEfectivo,1/scritico}

Here we see that   1/sfracEfect=1.437419940938726` is smaller than  1/(Subscript[S, ×])=1.6970562748477147` .  Therefore, the algorithm is unstable

Initial condition : u(x,t=0)= x(1-x)

uWA[nx_,0,alfa_]:=(nx+Nx)*deltax*(1-(nx+Nx)*deltax)

Table[{nx,uWA[nx,0,1]},{nx,-Nx,Nx}]//N

Boundary condition: u(0,t)=u(1,t)=0.

uL[t_]:=0;
uR[t_]:=0;

uWA[-Nx,m_,alfa_]:=uL[m*ht];
uWA[Nx,m_,alfa_]:=uR[m*ht];

The numerical algorithm  (see Eq.(19) of Ref[2])

RR[j_,m_,gamma_]:=Module[{kk},
kk=Sum[((1-lambda)*w[k+1,1-gamma]+lambda*w[k,1-gamma])*(uWA[j-1,m-1-k,gamma]-2*uWA[j,m-1-k,gamma]+uWA[j+1,m-1-k,gamma]),{k,0,m-1}];
uWA[j,m-1,gamma]+sfrac*kk]

uWA[j, m, gamma]  provides the value of the numeral solution at the point x_j, at time t_m for the case with anomalous diffusion exponent gamma

uWALista[m_,gamma_]:=uWALista[m,gamma]=LinearSolve[SparseArray[{Band[{1,2}]->Table[Which[j==Nx-1,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}],Band[{1,1}]->Table[Which[j==-Nx,1,j==Nx,1,True,(1+2*sfrac*w[0,1-gamma]*(1-lambda))],{j,-Nx,Nx}],Band[{2,1}]->Table[Which[j==-Nx,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}]},2*Nx+1],Table[Which[j==-Nx,uL[m*ht],j==Nx,uR[m*ht],True,RR[j,m,gamma]],{j,-Nx,Nx}]]

uWA[j_,m_,gamma_]:=uWA[j,m,gamma]=uWALista[m,gamma][[Nx+1+j]]

uWALista[2,gamma]

mfin=1;soluWA[mfin*ht,gamma]

mfin=10;soluWA[mfin*ht,gamma]

mfin=100;soluWA[mfin*ht,gamma]

mfin=100;ListPlot[soluWA[mfin*ht,gamma],Joined->True,PlotRange->All]

This bad result is due to the fact that the WA method is unstable for this case

Numerical solution of the FDE by means of the weighted-averaged method with weight factor λ=0.8  for  
	γ=1/2,   K=1,   ,   sfrac=0.55  ,   Δx=1/20 , Δt=1.89063×(10^(-6)) 
The algorithm is stable in this case

Here we clear the saved values of uWALista and uWA  (the evaluation of this cell is required if you are  going to work with several examples in a Mathematica session)

Clear[uWALista, uWA];
uWALista[m_,gamma_]:=uWALista[m,gamma]=LinearSolve[SparseArray[{Band[{1,2}]->Table[Which[j==Nx-1,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}],Band[{1,1}]->Table[Which[j==-Nx,1,j==Nx,1,True,(1+2*sfrac*w[0,1-gamma]*(1-lambda))],{j,-Nx,Nx}],Band[{2,1}]->Table[Which[j==-Nx,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}]},2*Nx+1],Table[Which[j==-Nx,uL[m*ht],j==Nx,uR[m*ht],True,RR[j,m,gamma]],{j,-Nx,Nx}]]
uWA[j_,m_,gamma_]:=uWA[j,m,gamma]=uWALista[m,gamma][[Nx+1+j]]

lambda=weight factor=λ

lambda=0.8;

2*Nx+1=number of spatial points

Nx=10;

Parameters of the example :  gamma = 1/2, Kg = 1, sfrac = 0.55, deltax = 1./(2*Nx)    
Notation:   gamma=γ,  ht=Δt  ,   Kg=K  ,   sfrac=(Subscript[S,  ]) ,  deltax=Δx

{gamma=0.5,deltax=1./(2*Nx),Kg=1.;sfrac=0.55,ht=(sfrac*deltax^2/Kg)^(1/gamma)}

According to Eq. (27) of Ref[2], the critical stability value for sfrac≡S=K(((Δt))^γ)/(((Δx))^2) is given by 
                          (Subscript[S, × ])= 1/(2*(2*lambda-1)*  w(-1,1-gamma)
As we use the BDF1 coefficients, we have
                                    w(-1,1-gamma)=2^(1-gamma)

The present numerical method is stable if   1/sfracEfect ≥  1/(Subscript[S, × ])with  sfracEfect=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2   
(see Eq. (25) of Ref[2])

{(Sin[(2*Nx-1)*Pi/(4*Nx)])^2//N;sfracEfectivo=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2;scritico=1/(2*(2*lambda-1)*2^(1-gamma));1/sfracEfectivo,1/scritico}

Here we see that   1/sfracEfect=1.8294435611947422` is larger than  1/(Subscript[S, ×])=1.6970562748477147` .  Therefore, the algorithm is stable

Initial condition : u(x,t=0)= x(1-x)

uWA[nx_,0,alfa_]:=(nx+Nx)*deltax*(1-(nx+Nx)*deltax)

Table[{nx,uWA[nx,0,1]},{nx,-Nx,Nx}]//N

Boundary condition: u(0,t)=u(1,t)=0.

uL[t_]:=0;
uR[t_]:=0;

uWA[-Nx,m_,alfa_]:=uL[m*ht];
uWA[Nx,m_,alfa_]:=uR[m*ht];

The numerical algorithm  (see Eq.(19) of Ref[2])

RR[j_,m_,gamma_]:=Module[{kk},
kk=Sum[((1-lambda)*w[k+1,1-gamma]+lambda*w[k,1-gamma])*(uWA[j-1,m-1-k,gamma]-2*uWA[j,m-1-k,gamma]+uWA[j+1,m-1-k,gamma]),{k,0,m-1}];
uWA[j,m-1,gamma]+sfrac*kk]

uWA[j, m, gamma]  provides the value of the numeral solution at the point x_j, at time t_m for the case with anomalous diffusion exponent gamma

uWALista[m_,gamma_]:=uWALista[m,gamma]=LinearSolve[SparseArray[{Band[{1,2}]->Table[Which[j==Nx-1,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}],Band[{1,1}]->Table[Which[j==-Nx,1,j==Nx,1,True,(1+2*sfrac*w[0,1-gamma]*(1-lambda))],{j,-Nx,Nx}],Band[{2,1}]->Table[Which[j==-Nx,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}]},2*Nx+1],Table[Which[j==-Nx,uL[m*ht],j==Nx,uR[m*ht],True,RR[j,m,gamma]],{j,-Nx,Nx}]]

uWA[j_,m_,gamma_]:=uWA[j,m,gamma]=uWALista[m,gamma][[Nx+1+j]]

uWALista[2,gamma]

mfin=1;soluWA[mfin*ht,gamma]

mfin=10;soluWA[mfin*ht,gamma]

mfin=100;soluWA[mfin*ht,gamma]

mfin=200;figNum=ListPlot[soluWA[mfin*ht,gamma],PlotRange->All,PlotStyle->PointSize[0.02]]

tv=200*ht;figExact=Plot[uExactx1mx[x,tv,1/2,1,20],{x,0,1},PlotStyle->{Red,Thick},PlotRange->All]

Show[figNum,figExact,PlotRange->All]

mfin=500;soluWA[mfin*ht,gamma]

mfin=1000;soluWA[mfin*ht,gamma]

mfin=1000;figNum=ListPlot[soluWA[mfin*ht,gamma],PlotRange->All,PlotStyle->PointSize[0.02]]

tv=1000*ht;figExact=Plot[uExactx1mx[x,tv,1/2,1,20],{x,0,1},PlotStyle->{Red,Thick},PlotRange->All]

Show[figNum,figExact,PlotRange->All]

Numerical solution of the FDE by means of the weighted-averaged method with weight factor λ=0, i.e.,  by means of a fully implicit method,  for  
	γ=1/2,   K=1,   ,   sfrac=0.7  ,   Δx=1/20 , Δt=6.25×(10^(-6)) .
Weighted-averaged methods with   λ≤1/2  are always stable;  see section 3.1 of  Ref[2]

Here we clear the saved values of uWALista and uWA  (the evaluation of this cell is required if you are  going to work with several examples in a Mathematica session)

Clear[uWALista, uWA];
uWALista[m_,gamma_]:=uWALista[m,gamma]=LinearSolve[SparseArray[{Band[{1,2}]->Table[Which[j==Nx-1,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}],Band[{1,1}]->Table[Which[j==-Nx,1,j==Nx,1,True,(1+2*sfrac*w[0,1-gamma]*(1-lambda))],{j,-Nx,Nx}],Band[{2,1}]->Table[Which[j==-Nx,0,True,-sfrac*w[0,1-gamma]*(1-lambda)],{j,-Nx,Nx-1}]},2*Nx+1],Table[Which[j==-Nx,uL[m*ht],j==Nx,uR[m*ht],True,RR[j,m,gamma]],{j,-Nx,Nx}]]
uWA[j_,m_,gamma_]:=uWA[j,m,gamma]=uWALista[m,gamma][[Nx+1+j]]

lambda=weight factor=λ

lambda=0;

2*Nx+1=number of spatial points

Nx=10;

Parameters of the example :  gamma = 1/2, Kg = 1, sfrac = 1, deltax = 1./(2*Nx)   -- -> always STABLE
Notation:   gamma=γ,  ht=Δt  ,   Kg=K  ,   sfrac=(Subscript[S,  ]) ,  deltax=Δx

{gamma=0.5,deltax=1./(2*Nx),Kg=1.;sfrac=1,ht=(sfrac*deltax^2/Kg)^(1/gamma)}

According to Eq. (27) of Ref[2], the critical stability value for sfrac≡S=K(((Δt))^γ)/(((Δx))^2) is given by 
                          (Subscript[S, × ])= 1/(2*(2*lambda-1)*  w(-1,1-gamma)
As we use the BDF1 coefficients, we have
                                    w(-1,1-gamma)=2^(1-gamma)

The present numerical method is stable if   1/sfracEfect ≥  1/(Subscript[S, × ])with  sfracEfect=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2   
(see Eq. (25) of Ref[2])

{(Sin[(2*Nx-1)*Pi/(4*Nx)])^2//N;sfracEfectivo=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2;scritico=1/(2*(2*lambda-1)*2^(1-gamma));1/sfracEfectivo,1/scritico}

Here we see that   1/sfracEfect=((Sec[(π/40)])^2) is larger than  1/(Subscript[S, ×])=-2.8284271247461903` .  Therefore, the algorithm is stable

Initial condition : u(x,t=0)= x(1-x)

uWA[nx_,0,alfa_]:=(nx+Nx)*deltax*(1-(nx+Nx)*deltax)

Table[{nx,uWA[nx,0,1]},{nx,-Nx,Nx}]//N

Boundary condition: u(0,t)=u(1,t)=0.

uL[t_]:=0;
uR[t_]:=0;

uWA[-Nx,m_,alfa_]:=uL[m*ht];
uWA[Nx,m_,alfa_]:=uR[m*ht];

uWALista[2,gamma]

mfin=1;soluWA[mfin*ht,gamma]

mfin=10;soluWA[mfin*ht,gamma]

mfin=100;soluWA[mfin*ht,gamma]

mfin=500;soluWA[mfin*ht,gamma]

mfin=500;figNum=ListPlot[soluWA[mfin*ht,gamma],PlotRange->All,PlotStyle->PointSize[0.02]]

tv=500*ht;figExact=Plot[uExactx1mx[x,tv,1/2,1,20],{x,0,1},PlotStyle->{Red,Thick},PlotRange->All]

Show[figNum,figExact,PlotRange->All]