m=2; A=Laplacian(m,2);             % setup Laplace equation as before
f=zeros(m*m,1); f(m:m:end)=1;
x=0:1/(m+1):1;
u=zeros(size(f)); U=zeros(m+2);
U(end,1:m+2)=1;
for k=1:4                          % do four Jacobi steps
  U(2:m+1,2:m+1)=reshape(u,m,m);  
  mesh(x,x,U);                     % plot current approximation
  xlabel('x');ylabel('y');
  u=u+(f-A*u)./diag(A);            % perform one Jacobi step
  pause
end