function root = secant ( f, a, b )
%
%  function root = secant ( f, a, b )
%
%  This is a version of the secant method.
%
%  F is the name of a MATLAB function file, of the form F(X).
%  A and B are two starting values at which to check F.
%
%  The program uses the secant method, starting from A and B,
%  and searches for a value at which the function F is near
%  enough to zero.  Up to 20 steps are taken.
%
FTOL = 0.00001;
XTOL = 0.00001;
NSTEP = 20;

fa = feval ( f, a ); 
if ( abs ( fa ) < FTOL )
  root = a;
  return
end

fb = feval ( f, b );
if ( abs ( fb ) < FTOL )
  root = b;
  return
end

for i = 1 : NSTEP

%
%  Can we take the step?
%
  if ( fb - fa == 0.0 ) 
    ' '
    'SECANT - Fatal error!'
    '  FA = FB, resulting in a zero denominator!'
    root = a;
    return
  end 
%
%  Take the step.
%
  c = ( fb * a - fa *  b ) / ( fb - fa );
  fc = feval ( f, c );
%
%  Is the new function value very small?
%
  if ( abs ( fc ) < FTOL )
    root = c;
    return
  end
%
%  Is the interval very small?
%
  if ( abs ( c - b ) < XTOL )
    root = 0.5 * ( b + c );
    return
  end
%
%  Update the data in preparation for the next step.
%
  a = b;
  b = c;

  fa = fb;
  fb = fc;

end
%
%  No convergence.
%
' '
'SECANT - Warning!'
'  Maximum number of steps, no convergence.'
'  Returning best estimate of root.'

root = 0.5 * ( a + b );