function [ x, it ] = newton ( f, fp, x, xtol, ftol, itmax )

%*****************************************************************************80
%
%% newton uses Newton's method to seek a root of a function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2019
%
%  Input:
%
%    function handle F: a pointer to a function.  This can be
%    an anonymous function or the name of a MATLAB M-file function.
%
%    function handle FP: a pointer to a function that evaluates the
%    derivative of F.
%    
%    real x: a starting point.
%
%    real xtol: the interval size tolerance.  
%    The test |x-xold| <= xtol is made.
%
%    real ftol: the residual tolerance.  
%    The test |f(x)| <= ftol is made.
%
%    integer itmax: the maximum number of steps allowed.
%    The test it <= itmax is made.
%
%  Output:
%
%    real x: the estimated root.
%
%    integer it: the number of steps.
%
  it = 0;
  new = 0.0;
  alpha = 0.0;

% fprintf ( 1, 'it = %d  x = %20.16g  f(x) = %g\n', it, x, f(x) );

  while ( true )

    it = it + 1;
 
    xold = x;
    x = x - f ( x ) / fp ( x );
%   fprintf ( 1, 'it = %d  x = %20.16g  f(x) = %g\n', it, x, f(x) );
    old = new;
    new = abs ( x - xold );

    alpha_old = alpha;
    if ( old ~= 0.0 )
      alpha = new / old;
      if ( 0 < alpha_old )
        r = log ( alpha ) / log ( alpha_old );
        fprintf ( 'it = %d: alpha = %g, log(alpha)= %g, r = %g\n', it, alpha, log(alpha), r );
      else
        fprintf ( 'it = %d: alpha = %g, log(alpha)= %g\n', it, alpha, log(alpha) );
      end
    end

    if ( abs ( x - xold ) <= xtol && abs ( f ( x ) ) <= ftol )
      return
    end

    if ( itmax <= it )
      return
    end
  end

  return
end

