nmax = 100;

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

tau = 1.1*tdiff;


t = 0:20*tau/100:50*tau;
x = 0:0.02:1;

for i = 1:length(t);
for j = 1:length(x);

u(j) = (tau/L)*(1.-exp(-t(i)/tau));

for n = 1:nmax;
term = (2/L)*cos(n*pi/2)/((n^2*pi^2*k/L^2) - (1./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
u(j) = u(j)/ tau
end
plot(x,u);
axis([0 1 0 5]);
h = gcf;
M(i) = getframe(h);
end

movie(M);

