% Heat diffusion in a wire

L = 1.0;
T = 0.3;
K = 1;
% m = 480;
% m=240;
 m=160
k = T/m;
n =20;
h = L/n;

u = zeros(n+1,m+1);

f = 0.0;
u0 = sin(pi*(0:n)*h)';
gLeft = 0.0;
gRight = 0.0;

u(:,1) = u0;

u(1,:) = gLeft;

u(n+1,:) = gRight;

x = (0:n)*h;
t = (0:m)*k;
alpha = K*k/h^2

for j=1:m
    for i=2:n
        u(i,j+1) = alpha*u(i-1,j) + (1 - 2*alpha)*u(i,j) + alpha*u(i+1,j) + f;
    end
end

%mesh(x,t,u')

for l=1:m+1
    plot(x,u(:,l),'b',x,exp(-pi^2*t(l))*sin(pi*x),'r');
    axis([0 1 0 1]);
    F(l) = getframe;
end

%movie(F)
