% m49_leapfrog2D.m - Leap frog for  u_tt = u_xx + u_yy
%      (cf. m47_BackwardEuler2D.m and m48_CrankNicolson2D.m)

% Grid and initial condition in the form of a Gaussian:
  J = input('grid parameter J (e.g. 10,30,60)? ');
  h = 1/J;
  s = (h:h:1-h)';
  k = .5*h;
  [xx,yy] = meshgrid(s,s);  
  x = xx(:); y = yy(:); 
  uold = exp(-50*((x-.5).^2+(y-.5).^2));
  u = uold;

% Matrices A and B:
  I = eye(J-1); II = eye((J-1)^2);         
  D = h^(-2)*toeplitz([-2 1 zeros(1,J-3)]);  
  L = kron(I,D) + kron(D,I);                 
  A = 2*II + k^2*L;

% Time-stepping:
  t = 0; shg
  while t<6
    uu = reshape(u,J-1,J-1);
    mesh(xx,yy,uu), colormap(1e-6*[1 1 1])
    title(sprintf('t = %4.2f',t), 'fontsize',20)
    axis([0 1 0 1 -1.2 1.2]), view(-20,30), drawnow
    t = t+k;
    unew = A*u - uold;
    uold = u; u = unew;
  end