function [u error i]= NewtonIteration(u0, AD, params, tol)
% Newton iteration to solve F(u)=0
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 
N=length(u0);
R=A*(1-u).*exp(-E./(u-p)); Fu=AD*u +R; Fu(1)=0; Fu(end)=0; % F(u), evaluated at the current state
for i=1:M
    % Jacobian matrix, evaluated at the current state
    Jr=-A*exp(-E./(u-p));  Jr=Jr+A*(1-u)*E.*exp(-E./(u-p))./(u-p).^2;
    Jr=spdiags(Jr, 0, N, N); % Jacobian of reactions
    J=AD+Jr;                 % Jacobian of system
    
    % Newton step
    v=J\(-Fu); 
    v(1)=0; v(end)=0;    % Boundary conditions to force v(0)=v(1)=0
    u=v+u;
    
    u(1)=1; u(end)=0;    % Boundary conditions to force u(0)=1 and u(1)=0
    u(u<0)=0; u(u>1)=1;  % force u to stay inside  0 <= u =< 1.
    
    R=A*(1-u).*exp(-E./(u-p)); Fu=AD*u +R; Fu(1)=0; Fu(end)=0;
    % When we are at the true solution, F(u)=0, 
    % and when we are near the true solution F(u) is small. 
    % So we use ||F(u)|| as an indication of the "error"
    % although strictly the error should be defined differently as u_approx-u_true.
    error(i)=norm(Fu);   
    if error(i)<tol
        break
    end
end