function u = poisson_jacobi(f)

% function u = poisson_jacobi(f)
% Jacobi iteration for - delta_h u = f on unit square

TOL = 1.e-6;
n = size(f,1);
ind = [2:n+1]; 

f = f / (n+1)^2;
u = zeros(n+2,n+2); v = u+2*TOL; 
while norm(v(:)-u(:),inf) > TOL*norm(u(:),inf)
   v = u; 
   u(ind,ind) = (f + v(ind-1,ind) + v(ind,ind-1) + ... 
                 v(ind+1,ind) + v(ind,ind+1) )/4;
end;

