function A = laplacian2(N)
% A = laplacian2(N)
% find SPARSE matrix A for finite element discretization of -Laplace u = f,
% using piecewise linear functions on a mesh with h=1/N.
%
% The matrix A has size M by M with M = (N-1)^2.

h = 1/N;       % meshsize
n = N-1;       % interior grid points are n by n grid
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
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)=4;  % A(k,k)   = a_diag
    if i>1
      m=m+1; iv(m)=k; jv(m)=k-1; av(m)=-1; % A(k,k-1) = a_left
    end
    if i<n
      m=m+1; iv(m)=k; jv(m)=k+1; av(m)=-1; % A(k,k+1) = a_right
    end
    if j>1
      m=m+1; iv(m)=k; jv(m)=k-n; av(m)=-1; % A(k,k-n) = a_down
    end
    if j<n
      m=m+1; iv(m)=k; jv(m)=k+n; av(m)=-1; % A(k,k+n) = a_up
    end
  end
end
iv = iv(1:m); jv = jv(1:m); av = av(1:m);

A = sparse(iv,jv,av,M,M);                   % assemble M by M sparse matrix from iv,jv,av


