Expand the domain for the 1D fractional heat equation

I have the Initial Value Problem for 1D fractional heat equation

with

And the exact solution being

So. I need to apply a numerical approximation to study it's rate of convergence

for

Here lies the problem. I have the program done. With the main issue being I can't expand the domain as asked without burning my computer

The following code takes 19 seconds to execute, if I expand the domain or reduce the step "h", the program takes too long to compute.

%format longE

tic
%Boundary conditions 
a = -100; %Need to be -5000
b = 100; %Need to be 5000
T = 1;

%Alpha value, in this program we only need alpha = 1
alpha = 1;

%Aproximation
h = .125; %Need to be 10^-2
tau = h^alpha/100;

%Boundary conditions 
f1 = @(x) 1./(1+x.^2);

%Exact solution
y = @(x,t) (t+1)./((t+1).^2 + x.^2);
 
%For the given values above, the result should be 0.0636. However, is the
%result I have is different

DFNLE(alpha,a,b,T,h,tau,f1,y);    


timeElapsed = toc

Main function

function DFNLE(alpha,a,b,T,h,tau,f1,y)

n = (b-a)/h -1;
n = floor(n);
x = a:h:b;
m = T / tau;
m = floor(m);
t = 0:tau:T;

w = zeros(n+2,1);
v1 = zeros(n,1);
v2 = zeros(n,1);

%Minimum jump
theta = h^(1/3);

%C constante
C = 2^alpha * gamma(1/2+alpha/2) / (sqrt(pi)*abs(gamma(-alpha/2)));

%Known values of u
u = zeros(n+2,m+1);
u(:,1) = f1(x);

%Loop to create the weights
 for k = 1:n+2
     if k*h > theta
        w(k) = C*h^(-alpha)/k^(alpha+1);
     end
 end
 
%CFL condition
if tau > 1/sum(w)
   fprintf('The methot does not converge for tau = %s', tau);
   return
end

%Main loop
for j = 1:m
    for i = 1:n+2

        for k = 1:n+2
 
%To calculate the sum, I separate it in two arrays, one that goes to the
%left and the other to the right
            if i+k <= n+2 
               v1(k) = (u(i+k,j)-u(i,j))*w(k);
            else
               v1(k) = -u(i,j)*w(k);
            
            end
                              
            if  i-k >= 1
               v2(k) = (u(i-k,j)-u(i,j))*w(k);   
            else 
               v2(k) = -u(i,j)*w(k);
            end
                           
        end  
        u(i,j+1) = u(i,j)+tau*(sum(v1)+sum(v2));  
    end
end

[X,t2] = meshgrid(x,t);

u = sparse(u);

%Error
err = max(max(abs(u' - y(X,t2))));
fprintf('Error for the explicit methord, epsilon = %d', err)

% Drawing 
% surf(X,t2,u');
% axis([-10 10 0 T 0 1])
% shading interp       
% title(['Explicit methor for alpha  = ', num2str(alpha)])
% xlabel('Variable x')
% ylabel('Variable t')
% zlabel('Función u') 

end