function [K,rhs,h] = laplace3d(n)

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

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

h = 1/(n+1);
K = K/(h*h);

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

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