% m47_BackwardEuler2D.m 
%       Backward Euler for  u_t = u_xx + u_yy on unit square
%       with periodic boundary data
%
%    xx, yy, uu represent data on a 2D mesh.
%
%    x, y, u represent the same data stretched to 1D vectors,
%            since that's what's needed to do linear algebra.

% Grid and initial condition in the form of a plus sign:
  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((abs(x-.5)<.4 & abs(y-.5)<.07) | (abs(x-.5)<.07 & abs(y-.5)<.4));

% The matrix B:
  I = speye(J); II = speye(J^2);               % identity matrices
  D = h^(-2)*toeplitz([-2 1 zeros(1,J-2)]);    % 1D Laplacian
 %D(1,J) = D(2,1); D(J,1) = D(1,2);            % periodic BC's
  L = kron(I,D) + kron(D,I);                   % 2D Laplacian
  B = II - k*L;
  spy(B), pause

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