% m56_waveeqcheb.m - Wave eq. u_tt = u_xx on [-1,1], Dirichlet BCs
%
%         t: leap frog.  x: Chebyshev spectral 

% Chebyshev grid and differentiation matrix:
  N = input('N? (e.g., 64 or 128) ');
  [D,x] = cheb(N);
  D2 = D^2;

% Time steps 0 and 1:
  k = 8/N^2;             % O(N^(-2)) time step needed for stability
  u = exp(-200*x.^2);
  hold off, shg
  plot(kron(x,[1 1]),1.95*[-1 1],'color',.55*[1 1 1])   % show grid
  hold on, plt = plot(x,u,'linewidth',3);
  axis([-1 1 -2 2])
  uold = u;
  u = exp(-200*(x-k).^2);

% Time-stepping:
  tmax = 8;
  nsteps = tmax/k;
  for n = 2:nsteps
    unew = 2*u - uold + k^2*D2*u;            
    uold = u; u = unew; u([1 end]) = 0;
 %  u(1) = -D(1,2:end)*u(2:end)/D(1,1);  % comment in for Neumann BC
    if rem(n,N^2/1024)==0
       set(plt,'ydata',u), drawnow
    end
  end