% demo of the heat equation u_t = u_{xx} with Dirichlet boundary conditions

 N = 8

% set up initial conditions
 x = chebfun('x',[0 1]);
 u0 = x.*(1-x);
 clear Psi

% find coefficients a_n(0)
 a0 = zeros(N,1);
 lambda = zeros(N,1);
 for n=1:N
    Psi{n} = sqrt(2)*sin(n*pi*x);
    lambda(n) = n^2*pi^2;
    a0(n) = (u0'*Psi{n})/(Psi{n}'*Psi{n}); 
 end 

% evolve solution in time
 
 tvec = [0:.0005:.005 .01:.01:.3];
 xx = linspace(0,1,100);
 U = zeros(length(xx),length(tvec));

 figure(1),  clf
 for j=1:length(tvec)
    t = tvec(j);
    u = 0*x;                                 % initialize u
    for n=1:N                                % sum first N terms in u(x,t)
       u = u + exp(-lambda(n)*t)*a0(n)*Psi{n}; 
    end
    plot(u,'b-','linewidth',2)
    U(:,j) = u(xx); 
    xlabel('$x$','fontsize',18,'interpreter','latex')
    ylabel('$u(x,t)$','fontsize',18,'interpreter','latex')
    title(sprintf('heat distribution at $t=%6.4f$',t),...
        'fontsize',18,'interpreter','latex')  
    ylim([0 .5])
    input('hit return')
 end

 figure(2), clf
 surf(xx, tvec,U')
 xlabel('$x$','interpreter','latex','fontsize',18)
 ylabel('$t$','interpreter','latex','fontsize',18)
 zlabel('$u(x,t)$','interpreter','latex','fontsize',18)

 figure(3), clf
 plt=waterfall(xx,tvec,U');
 set(plt,'edgecolor','k')
 xlabel('$x$','interpreter','latex','fontsize',18)
 ylabel('$t$','interpreter','latex','fontsize',18)
 zlabel('$u(x,t)$','interpreter','latex','fontsize',18)
