function [x, ithist, iflag] = fdnewtnd( f, x, tolf, tolx, maxit )
%
%  function [x, ithist, iflag] = fdnewtnd( f, x, tolf, tolx, maxit )
%
%  fdnewtnd attempts to compute a root of F using a finite difference 
%  Newton method
%
%  Input parameters:
%    f       name of a matlab function that evaluates 
%            f and its derivative.
%    x       initial iterate
%    tolf    stopping tolerance (optional. Default tolf = 1.e-7)
%            Newton's method stops if  ||F(x)||_2 < tolf
%    tolx    stopping tolerance (optional. Default tolx = 1.e-7)
%            Newton's method stops if  ||s||_2 < tolx, 
%            where s = -F'(x)^{-1} F(x)  is the Newton step.
%    maxit   maximum number of iterations (optional. Default maxit = 100)
%
%
%  Output parameters:
%    x       approximation of the solution. 
%    ithist  array with the iteration history
%            The i-th row of ithist contains  
%                        [it, norm(x), norm(F), norm(t*s), t]
%    ifag    return flag
%            iflag =  0  ||F(x)||_2 <= tolf 
%            iflag =  1  iteration terminated because maximum number of 
%                        iterations was reached. ||F(x)||_2 > tolf 
%%            iflag =  2  iteration terminated because no step size was found
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  March 11, 2004
%

% set tolerances if necessary
if( nargin <= 3 ) tolf = 1.e-7; tolx = 1.e-7; maxit = 100; end
if( nargin <= 4 ) tolx = 1.e-7; maxit = 100; end
if( nargin <= 5 ) maxit = 100; end


   
it       = 0;
iflag    = 0;
F        = feval(f, x);
% compute finite difference approximation of the Jacobian
h        = sqrt(eps)*norm(x);
Jac      = zeros(length(x),length(x));
for j = 1:length(x)
    xtmp = x;
    xtmp(j)  = xtmp(j) + h;
    Jac(:,j) = feval(f, xtmp);
    Jac(:,j) = (Jac(:,j) - F)/h;
end
    
normF    = norm(F);

s        = tolx*ones(size(x));
t        = 1;

while( it < maxit & norm(t*s) > tolx & normF > tolf )
   
   s     = - (Jac\F);
   
   t    = 1;
   xtmp = x+s;
   F    = feval(f, xtmp);
   while ( norm(F)^2 > (1-t*2.e-4)*normF^2 & norm(t*s) > tolx )
        t    = t/2;
        xtmp = x+t*s;
        F    = feval(f, xtmp);
   end

   ithist(it+1,:) = [it, norm(x), normF, norm(t*s), t];
   
   x  = xtmp;
   it = it+1;
   normF    = norm(F);
   % compute finite difference approximation of the Jacobian
   h        = sqrt(eps)*norm(x);
   Jac      = zeros(length(x),length(x));
   for j = 1:length(x)
       xtmp = x;
       xtmp(j)  = xtmp(j) + h;
       Jac(:,j) = feval(f, xtmp);
       Jac(:,j) = (Jac(:,j) - F)/h;
   end
   
end


% check why the FD-Newton method truncated and set iflag
if( norm(F) > tolf )
    % FD-Newton method truncated because the maximum number of iterations
    % was reached
    iflag = 1;
    return
elseif( norm(t*s) > tolx & t < 1 )
    iflag = 2;
    return
else
    % FD-Newton method truncated because norm(F) <= tolf
    % print info for last iteration
    ithist(it+1,:) = [it, norm(x), norm(F),0,0];
end