% demo of the wave equation u_{tt} = u_{xx}
% with Dirichlet boundary conditions, u(0,t)=u(1,t)=0
% and initial conditions u0 = x(1-x)*exp(2*x+.5*sin(12*pi*x)
% and initial velocity v0 = 0

% set up
 N = 20;
 x = chebfun('x',[0 1]);
 u0 = x.*(1-x).*exp(2*x+.5*sin(12*pi*x));   % initial position
 v0 = 0*x;
 clear Psi

% find coefficients a_n(0); construct eigenvalues and eigenfunctions
 b = zeros(N,1);
 d = 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;
    b(n) = (u0'*Psi{n})/(Psi{n}'*Psi{n});    
    c(n) = 0;
 end 

% evolve solution in time
 tvec = [0:.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)
       u = u + (b(n)*cos(sqrt(lambda(n))*t)+(c(n)/sqrt(lambda(n)))*sin(sqrt(lambda(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('wave $u(x,t)$ at $t=%6.4f$',t),...
        'fontsize',18,'interpreter','latex')  
    ylim([-1.5 1.5])
    input('hit return to advance time')
 end

% surface plot of solution
 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)

% waterfall plot of solution
 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)
