function [ W, A, B, IFLAG ] = rgfls ( F, A, B, XTOL, FTOL, NTOL )

%*****************************************************************************80
%
%% rgfls uses the regula falsi method to solve a nonlinear equation.
%
%***************  I N P U T  ****************
% F is the handle of the MATLAB function whose zero is sought
% A, B are the endpoints of the interval within the zero is sought
% XTOL desired length of the output interval
% FTOL desired size of F(W)
% NTOL upper bound on the number of iteration steps to be carried out
%**************  O U T P U T  ****************
% A, B the endpoints of the interval containing the zero
% W  best estimate of the zero
% IFLAG, an integer,
%      = -1, failure since F as the same size has the same sign at the 
%            input points A, B
%      = 0, termination because abs(A-B) <= XTOL
%      = 1, termination because abs(F(W) <= FTOL
%      = 2, termination after NTOL iteration steps 
%**************  M E T H O D  ****************
% The regula falsi algorithm 3.2 is used. This means that,
% at each step, linear interpolation between the points (A,FA) and
% (B,FB) is used, with FA*FB< 0, to get the new point (W,F(W))
% which replaces one of these in such a way that again FA*FB<0.
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
  FA = F(A);
  SIGNFA = sign(FA);
  FB = F(B); 
  if SIGNFA*FB>0
    fprintf('F is of the same sign at the two endpoints %e, %e\n',[A,B])
    IFLAG = -1;
    W = NaN;
    return
  end

  W = A;
  FW = FA;
  for n=1:NTOL
%
% check if interval is small enough
%
    if abs(A-B) <= XTOL
      IFLAG = 0;
      return
    end
%
% Check if F(W) is small enough
%
    if abs(FW) <= FTOL
      IFLAG = 1;
      return
    end
%
% Get new guess W by linear interpolation
%
    W = (FA*B-FB*A)/(FA-FB);
    PRVSFW = sign(FW);
    FW = F(W);
%
% Change to new interval
%
    if SIGNFA*FW > 0
      A = W;
      FA = FW; 
    else
      B = W;
      FB = FW;
    end

  end

  fprintf('no convergence in %d  iterations\n',NTOL)
  IFLAG = 2;

  return
end

