m=15;                             % resolution must be a power of 2 minus 1
l=5;                              % levels, l not bigger than L in m=2^L-1
h=1/(m+1);           
A=Laplacian(m,2)/h^2;             % construct the Laplacian on the fine grid
f=zeros(m*m,1);                   % right hand side for our model problem 
f(m:m:end)=1/h^2;
u=A\f;                            % exact solution
umg=zeros(size(f));               % multigrid initial guess 
U=zeros(m+2); U(end,1:m+2)=1;     % for plotting purposes
x=0:1/(m+1):1; y=x;               % mesh point vectors 
for n=0:10
  err(n+1)=max(max(abs(u-umg)));  % error in the infinity norm
  U(2:m+1,2:m+1)=reshape(umg,m,m);
  mesh(x,y,U)                     % plot current approximation
  xlabel('x');ylabel('y');
  pause                           % perform a V-cycle 
  umg=umg+VCycle(A,f-A*umg,zeros(size(f)),l);
end
