%Divide the grid
dx=0.1;dt=0.001;
r=0.1;a=1;
n=(1/dx)-1;
I=eye(n-1);
%Creat C matrix
for i=1:n-2
C(i,i+1)=1;
end
for i=2:n-1
C(i,i-1)=1;
end
%Accumulate D
D=(1-2*r*a)*I+r*a*C;
%Creat Vector
for i=2:n-1;
a(i)=D(i,i-1);
end
for i=1:n-1;
b(i)=D(i,i);
end
for i=1:n-2;
c(i)=D(i,i+1);
end
% Perform D=LU
u(1) = b(1);
for i = 2:n-1
l(i) = a(i)/u(i-1);
u(i) = b(i) - l(i) *c(i-1);
end
% Ly = d
y(1) = d(1);
for i = 2:n-1
d(i)=x(i-1);
y(i) = d(i)-l(i)*y(i-1);
end
% Ux = y
x(1) = y(1);
x(n-1) = y(n-1)/u(n-1);
for i = n-2:-1:1
x(i) = (y(i)-c(i)*x(i+1))/u(i); end