function result = brent ( f, xa, xb )
%
%  BRENT carries out Brent's method for seeking a real root of a nonlinear function.
%  This is a "stripped down" version with little error checking.
%
%  When the iteration gets going:
%    XB will be the latest iterate;
%    XA will be the previous value of XB;
%    XC will be a point with sign ( F(XC) ) = - sign ( F(XB) )
%
FATOL = 0.00001;
XATOL = 0.00001;
XRTOL = 0.00001;
ITMAX = 10;

it = 0;

fxa = feval ( f, xa );
fxb = feval ( f, xb );

xc = xa;
fxc = fxa;
d = xb - xa;
e = d;

format long

while ( it <= ITMAX )

  it = it + 1;

  [ xa, xb, xc; fxa, fxb, fxc ]

  if ( abs ( fxc ) < abs ( fxb ) ) 
    xa = xb;
    xb = xc;
    xc = xa;
    fxa = fxb;
    fxb = fxc;
    fxc = fxa;
  end

  xtol = 2.0 * XRTOL * abs ( xb ) + 0.5 * XATOL;

  xm = 0.5 * ( xc - xb );

  if ( abs ( xm ) <= xtol ) 
    'Interval small enough for convergence.'
    result = xb;
    return
  end

  if ( abs ( fxb ) <= FATOL )
    'Function small enough for convergence.'
    result = xb;
    return
  end
%
%  See if a bisection is forced.
%
  if ( abs ( e ) < xtol | abs ( fxa ) <= abs ( fxb ) )

    d = xm;
    e = d;

  else

    s = fxb / fxa;
%
%  Linear interpolation.
%
    if ( xa == xc )

      p = 2.0 * xm * s;
      q = 1.0 - s;
%
%  Inverse quadratic interpolation.
%
    else

      q = fxa / fxc;
      r = fxb / fxc;
      p = s * ( 2.0 * xm * q * ( q - r ) - ( xb - xa ) * ( r - 1.0 ) );
      q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );

    end

    if ( p > 0.0 )
      q = - q;
    else
      p = - p;
    end

    s = e;
    e = d;

    if ( 2.0 * p >= 3.0 * xm * q - abs ( xtol * q ) | p >= abs ( 0.5 * s * q ) )
      d = xm;
      e = d;
    else
      d = p / q;
    end

  end
%
%  Save old XB, FXB
%
  xa = xb;
  fxa = fxb;
%
%  Compute new XB, FXB,
%
  if ( abs ( d ) > xtol )
    xb = xb + d;
  elseif ( xm > 0.0 )
    xb = xb + xtol;
  elseif ( xm <= 0.0 )
    xb = xb - xtol;
  end

  fxb = feval ( f, xb );

  if ( sign ( fxb ) == sign ( fxc ) ) 
    xc = xa;
    fxc = fxa;
    d = xb - xa;
    e = d;
  end 

end

'Maximum number of steps taken.'
result = xc;