% demo of the wave equation u_{tt} = u_{xx} + sin(w*t) 
% with Dirichlet boundary conditions.
% Assume u(x,0) = u_t(x,0) = 0.
% Shows resonance when w is near an integer multiple of pi

% get the frequency at which to drive the wave equation
 w = input('enter w (near an integer multiple of pi): ');

% set up the eigenvalues and eigenfunctions
 N = 20;
 x = chebfun('x',[0 1]);
 clear Psi
 lambda = zeros(N,1);
 for n=1:N
    Psi{n} = sqrt(2)*sin(n*pi*x);
    lambda(n) = n^2*pi^2;
 end 

% evolve wave equation solution in time
 tvec = [0:.1:20];
 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)
       if abs(w-n*pi)<1e-14                  % w = n pi
          an = (sin(n*pi*t)-n*pi*t*cos(n*pi*t))/(2*n^2*pi^2);
       else
          an = (sin(w*t)-(w/(n*pi))*sin(n*pi*t))/(n^2*pi^2-w^2);
       end
       u = u + an*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([-2 2])
    pause(.1)
 end

% make a surface plot of the 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)

% make a waterfall plot of the 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)
