function [ A, B, IFLAG ] = bisect ( F, A, B, XTOL )

%*****************************************************************************80
%
%% bisect() shrinks the interval containing a zero of a function.
%
%  Discussion:
%
%    MATLAB version of SUBROUTINE BISECT on pages 81/2 of conte/deboor
%    returns the endpoints of the interval [A .. B] of length < XTOL
%    in the given interval [A .. B] that contains a zero of the function
%    whose function handle is F.
%
%    IFLAG is an integer
%       = -1, failure since F has the same sign at the given A,B
%       = 0,  abs(A-B)/2 <= XTOL
%       = 1,  abs(A-B)/2 is so small that its addition to (A+B)/2
%             makes no difference
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
  FA = F(A);

  if FA*F(B) > 0
    IFLAG = -1; 
    fprintf('F is of the same sign at the two endpoints %e, %e\n',[A,B]) 
    return
  end

  ERROR = abs(A-B);

  while ERROR > XTOL

    ERROR = ERROR/2;

    if ERROR <= XTOL
      IFLAG = 0;
      return
    end

    XM = (A+B)/2;
%
%  check for unreasonable error requirement
%
    if XM+ERROR==XM
      IFLAG = 1;
      return
    end

    FM = F(XM);
%
%  choose new interval
%
    if FA*FM >=0
      A = XM; 
      FA = FM;
    else
      B = XM;
    end

  end

  return
end
