function [ sol, it_hist, ierr ] = nsol ( x, f, tol, parms )

%*****************************************************************************80
%
%% nsol() is a locally convergent Newton solver for nonlinear systems.
%
%  Discussion:
%
%    This algorithm is a hybrid of Newton, Shamanskii, and Chord.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 April 2021
%
%  Author:
%
%    Tim Kelley
%
%  Reference:
%
%    Tim Kelley,
%    Iterative Methods for Linear and Nonlinear Equations,
%    SIAM, 2004,
%    ISBN: 0898713528,
%    LC: QA297.8.K45.
%
%  Input:
%
%    real X(N): the initial iterate.
%
%    function fx = F ( x ): MATLAB code that evaluates the function.
%
%    real TOL(2): [ atol, rtol ] error tolerances.
%    atol: absolute error tolerance.
%    rtol: relative error tolerance.
%
%    real PARMS(3): [ MAXIT, ISHAM, RSHAM ]
%    MAXIT: maximum number of iterations, default = 40.
%    ISHAM: the Jacobian matrix is recomputed and factored after ISHAM X updates.
%    default = 1000.
%    RSHAM: the Jacobian matrix is recomputed and factored if the ratio of
%    successive residual norms exceeds RSHAM, default = 0.5.
%		 isham = 1, rsham = 0 is Newton's method,
%		 isham = -1, rsham = 1 is the chord method,
%    isham = m, rsham = 1 is the Shamanskii method.
%
%  Output:
%
%    real SOL(N): the estimated solution.
%
%    real IT_HIST(num_iters): infinity norms of nonlinear residuals during iteration.
%
%    integer IERR: termination flag.
%    0: successful termination
%    1: after maxit iterations the termination criterion is not satsified,
%       or a ratio of successive residuals exceeded 1.
%
%  Local:
%
%    logical DEBUG: if TRUE, displays the iteration progress.
%
  debug = false;
  ierr = 0;
  maxit = 40;
  isham = 1000;
  rsham = 0.5;
  if ( nargin == 4 )
    maxit = parms(1);
    isham = parms(2);
    rsham = parms(3);
  end
  rtol = tol(2);
  atol = tol(1);
  it_hist = [];
  n = length ( x );
  fnrm = 1.0;
  itc = 0;
%
%  Evaluate f at the initial iterate
%
  f0 = feval ( f, x );
  fnrm = norm ( f0, inf );
  it_hist = [ it_hist, fnrm ];
  fnrmo = 1.0;
  itsham = isham;
%
%  Compute the stopping tolerance.
%
  stop_tol = atol + rtol * fnrm;
%
%  Iteration:
%
  while ( stop_tol < fnrm && itc < maxit )
%
%  Track the residual ratio.
%
    ratio = fnrm / fnrmo;
    outstat(itc+1, :) = [ itc, fnrm, ratio ];
    fnrmo = fnrm; 
    itc = itc + 1;
%
%  Update the Jacobian
%  * on the first iteration;
%  * every ISHAM iterates,
%  * if the ratio of successive residual norm is too large.
%
    if ( itc == 1 || rsham < ratio || itsham == 0 )
      itsham = isham;
      [ l, u ] = diffjac ( x, f, f0 );
    end
    itsham = itsham - 1;
%
%  Solve J * step = - f0
%
    tmp = - l \ f0;
    step = u \ tmp;

    xold = x;
    x = x + step;
    f0 = feval ( f, x );
    fnrm = norm ( f0, inf );
    it_hist = [ it_hist, fnrm ];
    ratio = fnrm / fnrmo;
    if ( debug )
      disp ( [ itc, fnrm, ratio ] )
    end
    outstat(itc+1,:) = [ itc, fnrm, ratio ];
%
%  If residual norms increase, terminate.
%
    if ( 1.0 <= ratio )
      ierr = 1;
      sol = xold;
      disp('nsol: increase in residual')
      disp(outstat)
      return;
    end

  end

  sol = x;
  if ( debug )
    disp(outstat)
  end
%
%  On failure, set the error flag
%
  if ( stop_tol < fnrm )
    ierr = 1;
  end

  return
end

