function v = mg5pt(f,u)
% MG5PT  multigrid iteration for 5 point Laplacian
%
% Usage: v = mg5pt(f,u)
%
%  f is a vector giving the right hand side
%  u is the initial guess, a grid function defined on the n x n interior grid points
%  output v is the improved solution after 1 V-cycle iteration
%
%  size of vectors must be n = 2^p-1 with p an integer
%
% To solve the 5 point Laplacian, we choose an initial guess for u (e.g., zero), and
% iterate several times: u=zeros(size(f)); for niter=1:niters; u=mg5pt(f,u); end
% Douglas N. Arnold, 2009-10-15

n = length(u);

if n == 1
  v = f/16;  % n=1 when h=1/2, so 4/h^2=16
else
  h = 1/(n+1);
  n2 = (n-1)/2;
% presmooth with Gauss-Seidel
  %  add zero values at boundary grid points
  U = [zeros(1,n+2);zeros(n,1),u,zeros(n,1);zeros(1,n+2)];
  for i=2:n+1
    for j=2:n+1
      U(i,j) = (U(i+1,j)+U(i-1,j)+U(i,j+1)+U(i,j-1)+h^2*f(i-1,j-1))/4;
    end
  end
  u = U(2:n+1,2:n+1);
% compute residual
  r = f -(4*u-U(1:n,2:n+1)-U(3:n+2,2:n+1)-U(2:n+1,1:n)-U(2:n+1,3:n+2))/h^2;
% transfer to coarser mesh
  r2 = ( 4*r(2:2:n-1,2:2:n-1) ...
        +2*(r(1:2:n-2,2:2:n-1)+r(3:2:n,2:2:n-1)+r(2:2:n-1,1:2:n-2)+r(2:2:n-1,3:2:n)) ...
        +r(1:2:n-2,1:2:n-2)+r(1:2:n-2,3:2:n)+r(3:2:n,1:2:n-2)+r(3:2:n,3:2:n) )/16;
% recursively apply multigrid iteration on coarse mesh
  e2 = mg5pt(r2,zeros(size(r2)));
% transfer solution to fine mesh
%  add boundary values
  E2 = [zeros(1,n2+2);zeros(n2,1),e2,zeros(n2,1);zeros(1,n2+2)];
  e = zeros(n,n);
  e(2:2:n-1,2:2:n-1) = e2;
  e(1:2:n,2:2:n-1) = (E2(1:n2+1,2:n2+1)+E2(2:n2+2,2:n2+1))/2;
  e(2:2:n-1,1:2:n) = (E2(2:n2+1,1:n2+1)+E2(2:n2+1,2:n2+2))/2;
  e(1:2:n,1:2:n) = (E2(1:n2+1,1:n2+1)+E2(1:n2+1,2:n2+2)+E2(2:n2+2,1:n2+1)+E2(2:n2+2,2:n2+2))/4;
% correct
  v = u + e;
% postsmooth with Gauss-Seidel in reverse order
  %  add zero valves at boundary grid points
  V = [zeros(1,n+2);zeros(n,1),v,zeros(n,1);zeros(1,n+2)];
  V(2:n+1,2:n+1) = v;
  for i=n+1:-1:2
    for j=n+1:-1:2
      V(i,j) = (V(i+1,j)+V(i-1,j)+V(i,j+1)+V(i,j-1)+h^2*f(i-1,j-1))/4;
    end
  end
  v = V(2:n+1,2:n+1);
end