f=zeros(m*m,1);                         % setup rhs of our model problem
f(m:m:end)=1;
u=A\f;                                  % exact solution
ud=zeros(size(f));                      % initial guess for ILU iterations
u0=ud;u001=ud;u002=ud;
errd=[]; err0=[]; err001=[]; err002=[]; % to store errors
uu=zeros(m+2); uu(m+2,1:m+2)=1;         % for the visualization
x=0:1/(m+1):1;
for n=0:10
  errd(n+1)=max(max(abs(u-ud)));          % compute errors
  err0(n+1)=max(max(abs(u-u0)));  
  err001(n+1)=max(max(abs(u-u001)));  
  err002(n+1)=max(max(abs(u-u002)));  
  uu(2:m+1,2:m+1)=reshape(ud,m,m);      % visualize iterates
%  subplot(2,2,1)
  mesh(x,x,uu,'LineWidth',2); title('LU') 
  xlabel('y');ylabel('x');
  set(gca,'FontSize',20)
  pause
  uu(2:m+1,2:m+1)=reshape(u0,m,m);
  subplot(2,2,2)
  mesh(x,x,uu,'LineWidth',2); title('ILU0') 
  xlabel('y');ylabel('x');
  set(gca,'FontSize',20)
  uu(2:m+1,2:m+1)=reshape(u001,m,m);
  subplot(2,2,3)
  mesh(x,x,uu,'LineWidth',2); title('ILUT(0.01)') 
  xlabel('y');ylabel('x');
  set(gca,'FontSize',20)
  uu(2:m+1,2:m+1)=reshape(u002,m,m);
  subplot(2,2,4)
  mesh(x,x,uu,'LineWidth',2); title('ILUT(0.002)') 
  xlabel('y');ylabel('x');
  set(gca,'FontSize',20)
  ud=ud+U\(L\(f-A*ud));                 % compute next iterate
  u0=u0+Ut0\(Lt0\(f-A*u0));
  u001=u001+Ut001\(Lt001\(f-A*u001));
  u002=u002+Ut002\(Lt002\(f-A*u002));
  pause
end
clf                                     % plot error curves
semilogy(0:10,errd,'-o',0:10,err0,'-s',0:10,err001,'-*',0:10,err002,'-x')
xlabel('Iteration'); ylabel('Error')
legend('LU','ILU0','ILUT(0.01)','ILUT(0.002)','Location','SouthEast')
set(gca,'FontSize',20)
axis([0 10 1e-15 1])