function [ x, x1, ierror, k ] = bisect ( x, x1, abserr, kmax, f )

%*****************************************************************************80
%
%% bisect() carries out the bisection method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, X1, two points defining the interval
%    in which the search will take place.  F(X) and F(X1) should have
%    opposite signs.
%
%    real ABSERR, an error tolerance.
%
%    integer KMAX, the maximum number of iterations allowed.
%
%    real external F, the name of the routine that
%    evaluates the function or its derivatives, of the form
%      value = function f ( x, ider )
%
%  Output:
%
%    real X, X1: X is the best estimate for the root, and X1 is
%    a recently computed point such that F changes sign in the interval
%    between X and X1.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%

%
%  Initialization.
%
  ierror = 0;
  k = 0;
%
%  Evaluate the function at the starting points.
%
  fx = f ( x, 0 );
  fx1 = f ( x1, 0 );
%
%  Set XPOS and XNEG to the X values for which F(X) is positive
%  and negative, respectively.
%
  if ( 0.0 <= fx && fx1 <= 0.0 )

  elseif ( fx <= 0.0 && 0.0 <= fx1 )
    t = x;
    x = x1;
    x1 = t;
    t = fx;
    fx = fx1;
    fx1 = t;
  else
    ierror = 1;
    return
  end
%
%  Iteration loop:
%
  while ( true )
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      break;
    end

    k = k + 1;

    if ( kmax < k )
      ierror = 2;
      break;
    end
%
%  Set the increment.
%
    dx = 0.5 * ( x1 - x );
%
%  Update the iterate and function values.
%
    x2 = x + dx;
    fx2 = f ( x2, 0 );

    if ( 0.0 <= fx2 );
      x = x2;
      fx = fx2;
    else
      x1 = x2;
      fx1 = fx2;
    end

  end

  return
end
