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

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

% 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
    U = fft(u);
    W = -[0:N/2 -N/2+1:-1]'.^2.*U;    % 2nd-order diff. via FFT
    w = real(ifft(W));
    unew = 2*u - uold + k^2*w;            
    uold = u; u = unew;
    if rem(n,N/8)==0
       set(plt,'ydata',u), drawnow
    end
  end