% m48_CrankNicolson2D.m
%   Crank-Nicolson for  u_t = u_xx + u_yy
%   (cf. m47_BackwardEuler2D.m)

% Grid and initial condition in the form of a circle:
  J = input('grid parameter J (e.g. 10,30,90)? ');
  h = 1/J;
  s = (h:h:1)';
  k = .0002;
  [xx,yy] = meshgrid(s,s);  
  x = xx(:); y = yy(:); 
  u = double(sqrt((x-.5).^2+(y-.5).^2) < .25);

% Matrices A and B:
  I = speye(J); II = speye(J^2);         
  D = h^(-2)*toeplitz([-2 1 zeros(1,J-2)]);  
 %D(1,J) = D(2,1); D(J,1) = D(1,2);         % periodic BC's
  L = kron(I,D) + kron(D,I);                 
  A = II + k*L/2;
  B = II - k*L/2;

% Time-stepping:
  t = 0;
  while 1
    uu = reshape(u,J,J);
    surfl(xx,yy,uu), colormap(hot)
    title(sprintf('t = %6.4f',t), 'fontsize',20)
    axis([0 1 0 1 0 2]), view(-20,30), pause
    t = t+k;
    u = B\(A*u);      % removing the parens is catastrophic!
  end