function A = spgridmat(n,aD,aW,aE,aS,aN)
% A = spgridmat(n,aD,aW,aE,aS,aN)
% find matrix A for problem on n by n grid:
% consider n^2 unknowns u_1,...,u_{n^2} on n by n grid
% arranged like this:
%                     ...   u_{n^2}         
%     :       :             :
%     :       :             :
%    u_{n+1}  u_{n+2} ...   u_{2n} 
%    u_1      u_2     ...   u_n
%
% each unknown u has up to four neighbors uW, uE, uS, uN in West, East, South, North
% for each unknown u we have an equation 
%   aD*u + aW*uW + aE*uE + aS*uS + aN*uN = ...
% (at boundary some terms may be missing)

M = n^2;       % number of interior grid points. Matrix A is of size M by M
nzmax = 5*M;   % number of nonzero matrix elements is <= 5*M (at most 5 nonzero elements per row)
iv = zeros(nzmax,1); jv = zeros(nzmax,1); av = zeros(nzmax,1); % initialize vectors
% describe nonzero matrix elements by vectors iv, jv, av of same length
% iv contains i values, jv contains j values, av contains matrix entries
m = 0;
for i=1:n      % for columns 1 to n
  for j=1:n    %   for rows 1 to n
    k = i + n*(j-1);      % k is number of point in row j, column i
    m=m+1;   iv(m)=k; jv(m)=k;   av(m)=aD;   % diagonal entry
    if i>1
      m=m+1; iv(m)=k; jv(m)=k-1; av(m)=aW;   % k-1 is number of point to the West
    end
    if i<n
      m=m+1; iv(m)=k; jv(m)=k+1; av(m)=aE;   % k+1 is number of point to the East
    end
    if j>1
      m=m+1; iv(m)=k; jv(m)=k-n; av(m)=aS;   % k-n is number of point to the South
    end
    if j<n
      m=m+1; iv(m)=k; jv(m)=k+n; av(m)=aN;   % k+n is number of point to the North
    end
  end
end
A = sparse(iv(1:m),jv(1:m),av(1:m)); % assemble sparse matrix from iv,jv,av

