function [u error i]= FPIteration(u0, AD, params, tol)
% Fixed point iteration 
AD(1,1)=1; AD(end,end)=1; % boundary conditions 
A=params.A; E=params.E; p=params.p; 
M=35;     % maximum iterations
u=u0;     % initial guess 
for i=1:M
    R=A*(1-u).*exp(-E./(u-p));
    R(1)=-1; R(end)=0; % boundary conditions to force u(0)=1 and u(1)=0
    u=AD\(-R);
    
    % evaluate right hand side at new guess
    R=A*(1-u).*exp(-E./(u-p));
    Fu=AD*u +R; Fu(1)=0; Fu(end)=0;
    u(u<0)=0; u(u>1)=1;  % force u to stay inside  0 <= u =< 1.
    error(i)=norm(Fu);   
    if error(i)<tol
        break
    end
end