function [K,rhs,h,A] = laplace2d(n)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [K,rhs,h] = laplace2d(n)
% Forms matrix K for finite difference approximation K u = rhs
% to Laplacian - u_{xx} - u_{yy} on an n x n grid
% and generates a smooth right-hand side rhs.
% The output parameter h is the mesh length: h = 1/(n+1)
% and A is the matrix for the 1-d problem.
%
% Dianne O'Leary 04-2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h = 1/(n+1);
A = spdiags([-ones(n,1),2*ones(n,1),-ones(n,1)],-1:1,n,n)/(h*h);
en = speye(n);
K = kron(A,en) + kron(en,A);

for i=1:n,
  x = i*h;
  for j=1:n,
    y = j*h;
    f(i,j) = (x^2+y^2)*cos(x*y*pi);
  end
end

rhs = reshape(f,n^2,1);