function A = setupA (N);
%setupA system matrix for finite difference discretization on the unit square
%   setupA sets up the system matrix A for finite difference discretization
%   of the Poisson equation on the unit square using the classical
%   five-point stencil for the Laplace operator with N+1 subintervals
%   in each direction. A is a N-by-N block tridiagonal matrix of
%   N-by-N block matrices. There are -1's in the off-diagonal entries and
%   4's on the diagonal. Note that the factor 1/h^2 is *not* included in A.

%   2001-2002 by Matthias K. Gobbert

  % ... auxiliary variables:
  n = N^2;
  E = ones (N, N);

  % ... set up diagonals of A in columns of matrix C using `row orientation':
  C = zeros (n, 5);
  V = E; V( : , 1 ) = 0; C(:,1) = -1.0 * V(:);
  V = E; V( 1 , : ) = 0; C(:,2) = -1.0 * V(:);
  V = E;                 C(:,3) = +4.0 * V(:);
  V = E; V(end, : ) = 0; C(:,4) = -1.0 * V(:);
  V = E; V( : ,end) = 0; C(:,5) = -1.0 * V(:);

  % ... set positions of diagonals in A and convert C to `column orientation':
  d = [-N -1 0 +1 +N];
  for i = 1 : length(d)
    I = [ max(1,1-d(i)) : min(n,n-d(i)) ];
    J = [ max(1,1+d(i)) : min(n,n+d(i)) ];
    B(J,i) = C(I,i);
  end;

  % ... set up A in sparse matrix storage:
  A = spdiags (B, d, n, n);

