% m51_waveeqFourier.m - Wave eq. u_tt = u_xx on [-pi,pi], periodic BCs
%
%         t: leap frog.  x: spectral, via matrices

% Grid:
  N = input('N? (e.g., 64 or 256) ');
  h = 2*pi/N;
  x = -pi+h*(1:N)';
  k = .1*h;          

% Spectral differentiation matrix:
% (see p. 23 of Trefethen, Spectral Methods in Matlab, SIAM 2000)
  c = [-pi^2/(3*h^2)-1/6 -.5*(-1).^(1:N-1)./sin(h*(1:N-1)/2).^2];
  D = toeplitz(c);

% Time steps 0 and 1:
  u = exp(-20*x.^2);
  hold off, subplot(2,1,1), shg
  plt = plot(x,u,'linewidth',4);
  axis([-pi pi -1.2 1.2]), grid
  uold = u;
  u = exp(-20*(x-k).^2);

% Time-stepping:
  tmax = 8*pi;
  nsteps = tmax/k;
  for n = 2:nsteps
    unew = 2*u - uold + k^2*D*u;   
    uold = u; u = unew;
    if rem(n,N/8)==0
       set(plt,'ydata',u), drawnow
    end
  end