nmax = 1000;

L = 1;
k = 1.0;
tdiff = L^2/(pi^2*k);

tau = tdiff/16.1


t = 0:0.001:1;
x = 0:0.02:1;



for i = 1:length(t);
for j = 1:length(x);
u(j) = 0.;
for n = 0:nmax;

term = cos(n*pi/2)/((n^2*pi^2*k/L^2) - (1./tau))/tau;
u(j) = u(j) + (exp(-t(i)/tau) - exp(-n^2*pi^2*k*t(i)/L^2) )*cos(n*pi*x(j)/L)*term;
end
end
plot(x,u);
axis([0 1 0 5]);
h = gcf;
M(i) = getframe(h);
end

%movie(M);

