function [ ido, x ] = newton_rc ( ido, x, fx )

%*****************************************************************************80
%
%% newton_rc() solves a small system of nonlinear equations.
%
%  Discussion:
%
%    NEWTON_RC uses Newton's method, with the jacobian approximated via
%    finite differences.
%
%    NEWTON_RC uses "reverse communication".  That is, NEWTON_RC does not call 
%    any user subroutines.  Instead, the user repeatedly calls NEWTON_RC, 
%    and the user and NEWTON_RC communicate via values of the parameter IDO.
%
%    To begin using NEWTON_RC, set IDO=0, set X to an approximate root,
%    set FX(1) to a small positive tolerance, and call NEWTON_RC.
%
%    NEWTON_RC will return with IDO=2.  Evaluate the residual at the
%    point X that is returned by NEWTON_RC, store that value in FX,
%    and call NEWTON_RC back with IDO=2.
%
%    NEWTON_RC will return with IDO=1.  Again, evaluate the residual
%    at the point X returned by NEWTON_RC, store that value in FX<
%    and call NEWTON_RC back with IDO=1.
%
%    This process will be repeated until NEWTON_RC returns with
%    IDO=0, meaning that the value X is a good approximation to
%    the root, or IDO=-1, in which case the algorithm failed.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    07 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer IDO, indicates the user's action.
%    * 0, this is a new problem, with N, X, and FX supplied.
%    * 1 or 2, FX evaluated at the X returned by previous call.
%
%    integer N, the number of equations.
%
%    real X(N), the current iterate.
%
%    real FX(N): the function value at the current X.
%
%  Output:
%
%    integer IDO.
%    * 0, the output value X is a good approximation to the solution.
%    * 1, please evaluate F(X).  (Jacobian is being approximated).
%    * 2, please evaluate F(X).
%    * -1, the algorithm failed.
%
%   real X(N), the current point being considered.
%
  persistent delxj
  persistent epsilon
  persistent fprime
  persistent fxcall
  persistent fxnrm
  persistent fxnrm0
  persistent fxold
  persistent j
  persistent ncall
  persistent xold

  x = x(:);
  fx = fx(:);
  n = size ( x, 1 );
%
%  IDO=0, user is starting a new problem.
%
  if ( ido == 0 )

    delxj = 0.0;
    epsilon = sqrt ( sqrt ( eps ) );
    fprime = zeros ( n, n );
    fxcall = 0;
    fxnrm = norm ( fx );
    fxnrm0 = fxnrm;
    fxold = fx;
    j = 0;
    ncall = 0;
    xold = x;

    ido = 2;
%
%  IDO=1, user is returning F(X+delX(J)) for a jacobian estimation
%
  elseif ( ido == 1 )

    if ( 1000 < ncall )
      fprintf ( 1, 'newton_rc: too many calls!\n' );
      ido = -1;
      return
    end

    ncall = ncall + 1;
    fprime(1:n,j) = ( fx(1:n) - fxold(1:n) ) / delxj;

    x(j) = xold;

    if ( j < n )

      j = j + 1;
      delxj = epsilon * ( abs ( x(j) ) + 1.0 );
      xold = x(j);
      x(j) = x(j) + delxj;

    else

      r = rcond ( fprime );
      if ( r < 1.0e-15 )
        fprintf ( 1, 'newton_rc: jacobian is ill-conditioned.\n' );
      end

      delx = - fprime \ fx;
      x(1:n) = x(1:n) + delx(1:n);

      ido = 2;

    end
%
%  IDO = 2, user is returning F(X) for a convergence evaluation.
%
  elseif ( ido == 2 )

    if ( 1000 < ncall )
      fprintf ( 1, 'newton_rc: too many calls!\n' );
      ido = -1;
      return
    end

    ncall = ncall + 1;
    fxcall = fxcall + 1;

    fxnrm_previous = fxnrm;

    xnrm = norm ( x );
    fxnrm = norm ( fx );

    if ( fxnrm <= eps * ( fxnrm0 + 1.0 ) )
      ido = 0;
      return
    end

    if ( 15 < fxcall )
      if ( 0.95 * fxnrm_previous < fxnrm )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'newton_rc - Warning\n' );
        fprintf ( 1, '  Previous |FX| = %g\n', fxnrm_previous );
        fprintf ( 1, '  Latest   |FX| = %g\n', fxnrm );
        fprintf ( 1, '  Convergence has slowed down.\n' );
        ido = -1;
        return
      end
    end
%
%  Prepare to move to jacobian estimation.
%
    ido = 1;
    fxold(1:n) = fx(1:n);
    j = 1;
    delxj = epsilon * ( abs ( x(j) ) + 1.0 );
    xold = x(j);
    x(j) = x(j) + delxj;
%
%  Unexpected value of IDO.
%
  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'newton_rc(): Fatal error.\n' );
    fprintf ( 1, '  Unexpected value of IDO = %d\n', ido );
    error ( 'newton_rc(): Fatal error.' )

  end

  return
end

