% demo of the heat equation u_t = u_{xx} + f(x) with Dirichlet boundary conditions

% set up
 N = 50;
 x = chebfun('x',[0 1]);
 u0 = x.*(1-x).*exp(2*x+.5*sin(12*pi*x));   % initial condition
 f  = 100*x.^1.*sin(10*pi*x);               % forcing function (t independent)
 clear Psi

% find coefficients a_n(0)
 a0 = zeros(N,1);
 lambda = zeros(N,1);
 ulim = 0*x;
 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}); 
    c(n) = (f'*Psi{n})/(Psi{n}'*Psi{n}); 
    ulim = ulim + (c(n)/lambda(n))*Psi{n};
 end 
 

% evolve solution in time
 
 tvec = [0:.0005:.005 .01:.01:.5];
 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)
       ant = exp(-lambda(n)*t)*a0(n)+((1-exp(-lambda(n)*t))/lambda(n))*c(n);
       u = u + ant*Psi{n}; 
    end
    clf
    plot(ulim,'r-','linewidth',2)
    hold on
    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([-.25 1.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)

