close all; clear; clc;
alpha = 1
N  = 100;
Lx = 1;
dx = Lx/(N-1);
x  = 0:dx:Lx;
M = 1000;
tf = 1;
dt = tf/(M-1);
t = 0:dt:tf;
lambda = alpha*dt/(dx^2);
f  = @(x) sin(pi*x/Lx);
fanal = @(x,t) sin(pi*x/Lx)*(exp(1))^(-[pi^2]*alpha*t/Lx);

u = zeros(N,M);
u_anal = zeros(N,M);

u(:,1) =  f(x);
u_anal(:,1)= fanal(x,0);

for k= 1:M-1
     for i= 2:N-1
       u(i,k+1) = (1-2*lambda)*u(i,k) + lambda*[(u(i-1,k))+ u(i+1,k)];
       u_anal(i,k+1)=fanal(x(i),t(k));
     end
end


for i = 1:6
    figure
    plot(x,u_anal(:,i),'-r',x,u(:,i),'--b');
    a = ylabel('Temperature');
    set(a,'Fontsize',14);
    a = xlabel('x');
    set(a,'Fontsize',14);
    a=title(['Temperature distribution at t = ' num2str(i-1)])
    legend('Analytic distribution','Estimated distribution')
    set(a,'Fontsize',16);
    xlim([0 1]);
    grid;
end