% m50_waveeq.m - Wave eq. u_tt = u_xx on [-pi,pi], periodic BCs
%
%         t: leap frog   x: finite diffs of orders 2,4 or 6

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

% Sparse finite differencing matrices of orders 2, 4 or 6:
% (See Fornberg, Math. Comp. 51 (1988), 699-706.)
  I = speye(N);
  h2 = h^(-2);
  shift = [N 1:N-1];
  I1 = I(:,shift); I2 = I1(:,shift); I3 = I2(:,shift);
  I1 = I1+I1'; I2 = I2+I2'; I3 = I3+I3';
  order = input('order 2, 4, or 6? ');
  switch order
    case 2, D = h2*( I1 - 2*I );
    case 4, D = h2*( -(1/12)*I2 + (4/3)*I1 - (5/2)*I );
    case 6, D = h2*( (1/90)*I3 - (3/20)*I2 + (3/2)*I1 -(49/18)*I);
  end

% Time steps 0 and 1:
  u = exp(-20*x.^2);
  hold off, subplot(2,1,1), shg
  plt = plot(x,u,'linewidth',2);
  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