BTCS for 1D heat problem

Hello all

I am working on the probelem of solving the heat equation with 3 methods of exact,FTCS, and BTCS. The fisrt two methods work fine. But the result for BTCS (solving uisng thomas algorithm) is not correct. Can anyone help with troubleshooting this?

    
clear;clc;;close all;
%setting up the parameters
Nmax=500; %max of Fourier modes for fine case
U0=2;
tf=10;
dt=0.5;
dx=0.1; 
Kappa=5*10^-3;
x = (-3:dx:3);
n=length (x);
alpha=Kappa*dt./dx.^2;
% problem initialization
    for j=1:n
  if (1<x(j)) || (x(j)<-1)
    Tinit(j)=0;
  else Tinit(j)=U0;
  end
 end
   %********************************************************************************
 
     %Calculating F exact
 for k=-3:dx:3
     fexact=0.5*U0*(erf((1-x)/(2*sqrt(Kappa*tf)))-erf(-(x+1)/(2*sqrt(Kappa*tf))));
  end
     %************************************************************************************
   %Calculating T_ftcs
   T0=Tinit;
 for l=1:dt:tf
for i=2:n-1
T_ftcs(i) = T0(i) + alpha*(T0(i+1)-2*T0(i)+T0(i-1));
T_ftcs(1)=T0(1) + alpha*(T0(2)-2*T0(1)+0);
T_ftcs(n) = T0(n) + alpha*(0-2*T0(n)+T0(n-1));
end 
T0=T_ftcs;
 end
 
%%****************************************************************************
 % calculating f_BTCS
 % Creat matrix A
e= ones(n-2,1);
Ac =spdiags([-alpha*e (1+2*alpha)*e -alpha*e],-1:1,n-2,n-2);
%initials
T_btcs=Tinit;
%start loop in time
 for s=1:dt:tf
b =T_btcs(2:n-1);
b(1)= b(1)+alpha*T_btcs(1);
b(end)= b(end)+alpha*T_btcs(end);
p=length(b);
%solving linear system with Thomas algorithm

  for k = 1:p-1
i = k +1;
l(i,k)=Ac(i,k)/Ac(k,k);
for j = k :k+1
Ac(i,j) = Ac(i,j)- l(i,k)*Ac(k,j);
end
b(i) = b(i)-l(i,k)*b(k);
  end
   %apply backward substitution
for k = p:-1:1
x(k)=b(k);
for j=k+1: min(p,k+1)
    x(k)=x(k)-Ac(k,j)*x(j);
end
x(k)=x(k)/Ac(k,k);

end
T_btcs(k)=x(k);
 end
 
 
 figure ('units', 'normalized','position',[0.55 0.1 .45 .45])
plot(x,fexact,'o',x,T_ftcs,'-d',x,T_btcs,'--')
set(gca,'fontsize',26)
ylim([0,3])
xlabel ('x')