function [x,info] = newton(fct,x,par)

% function [x,info] = newton(fct,x,par)
% solution of a nonlinear system f(x) = 0 
% with damped Newton iteration 
% 
% fct: name of function [f,df] = fct(x)
%      evaluating f and the Jacobi matrix df
% x:   start vector and solution
% par.tol: tolerance
% par.nit: maximal number of iterations
% par.step: minimal damping parameter
% info.msg: error message 
% info.it: number of iterations
% info.err: norm of f(x)

% defaults 
if ~exist('par'); 
   par = struct('tol',[],'nit',[],'step',[]);
end;
if isempty(par.tol);
   par.tol = 10^(-10);
end;
if isempty(par.nit);
   par.nit = 50; 
end;
if isempty(par.step);
   par.step = 1/512; 
end;

% initialization
s = 1; 

for k=1:par.nit;
   [f,df] = feval(fct,x);

   % test for convergence
   info.err = norm(f,inf);
   if info.err <= par.tol; 
      info.it = k; 
      return; 
   end;

   % computation of increment
   [Q,R] = qr(df);
   if min(abs(diag(R))) < par.tol*norm(R,inf); 
      info.msg = 'near singular Jacobi matrix'; 
      return;
   end;
   dx = rsolve(R,Q'*f);
   
   % checking increment and damping
   flag_s = 0;
   while s > par.step;
      y = x - s*dx;
      g = feval(fct,y);
      dy = rsolve(R,Q'*g);
      if norm(dy,inf) <= (1-s/2) * norm(dx,inf);
	 x = y;
	 if flag_s == 0; 
	    s = min(1,2*s);
         end;
	 break; 
      else; 
         s = s/2;
	 flag_s = 1;
      end;
   end; 
   if s <= par.step; 
      info.msg = 'damping parameter too small';
      return;
   end;
end;
info.msg = 'maximal number of iterations exceeded';


% local function 

function d = rsolve(R,f) 
% solves upper triangular system 

n = length(f);
d = [zeros(n-1,1); f(n)/R(n,n)];
for k=n-1:-1:1;
   d(k) = (f(k) - R(k,k+1:n)*d(k+1:n))/R(k,k);
end; 

