function x = secant(f,x0,x1)
% SECANT   Secant method for a scalar equation.
% Input:
%   f        function that outputs values of f(x)
%   x0,x1    initial root approximations
% Output       
%   x        vector of root approximations (last is best)

% Operating parameters.
funtol = 100*eps;  xtol = 100*eps;  maxiter = 40;

x = [x0 x1];
dx = Inf;  fn0 = f(x0);  fn1 = f(x1);
n = 2;

while (abs(dx) > xtol) & (abs(fn1) > funtol)
  dx = -fn1 * (x(n)-x(n-1)) / (fn1-fn0);   % secant step
  x(n+1) = x(n) + dx;

  n = n+1;
  if n==maxiter
    warning('Maximum number of iterations reached.')
    break
  end

  fn0 = fn1;
  fn1 = f(x(n));
end