m=15;                                   % number of gridpoints  
A=Laplacian(m,2);                       % five point Laplacian
f=zeros(m*m,1); f(m:m:end)=1;           % put bc into the rhs
u=A\f;                                  % solve by sparse Gaussian elimination
h=1/(m+1); x=0:h:1; y=x;                % mesh point vectors
un=zeros(size(f));                      % initial guess
UU=zeros(m+2); UU(end,1:m+2)=1;         % for plotting
[L,U]=ILU0(A);
for n=0:10
  err(n+1)=max(max(abs(u-un)));         % compute error 
  UU(2:m+1,2:m+1)=reshape(un,m,m);
  mesh(x,y,UU)
  xlabel('x');ylabel('y');
  pause
  un=un+U\(L\(f-A*un));
end
