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

%*****************************************************************************80
%
%% newton_system2 uses Newton's method to solve a system of nonlinear equations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 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
%    jacobian matrix of F.
%    
%    real x(n): 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(n): the estimated root.
%
%    integer it: the number of steps.
%
  global print_f_norm
  global print_r
  global print_x

  it = 0;
  dx = 0.0;
  alpha = 0.0;

  if ( print_x )
    fprintf ( 1, 'it = %d, x = (', it )
    for i = 1 : length ( x )
      fprintf ( 1, ' %20.16g', x(i) );
    end
    fprintf ( 1, ' )\n' );
  end

  if ( print_f_norm )
    fprintf ( 1, 'it = %d, ||f(x)|| = %20.6f\n', it, norm ( f(x) ) )
  end

  while ( true )

    it = it + 1;
 
    dx_old = dx;
    dx = - fp ( x ) \ f ( x );
    x = x + dx;

    if ( print_x )
      fprintf ( 1, 'it = %d, x = (', it )
      for i = 1 : length ( x )
        fprintf ( 1, ' %20.16g', x(i) );
      end
      fprintf ( 1, ' )\n' );
    end

    if ( print_f_norm )
      fprintf ( 1, 'it = %d, ||f(x)|| = %20.6f\n', it, norm ( f(x) ) )
    end

    alpha_old = alpha;
    if ( norm ( dx_old ) ~= 0.0 )
      alpha = norm ( dx ) / norm ( dx_old );
      if ( 0 < alpha_old )
        r = log ( alpha ) / log ( alpha_old );
        if ( print_r )
          fprintf ( 'it = %d: alpha = %g, log(alpha)= %g, r = %g\n', ...
            it, alpha, log(alpha), r );
        end
      else
        if ( print_r )
          fprintf ( 'it = %d: alpha = %g, log(alpha)= %g\n', ...
            it, alpha, log(alpha) );
        end
      end
    end

    if ( norm ( dx ) <= xtol && norm ( f ( x ) ) <= ftol )
      return
    end

    if ( itmax <= it )
      return
    end
  end

  return
end

