function x = bisect(f,a,b)
% BISECT   Bisection method for a scalar equation.
% Input:
%   f        function that outputs values of f(x)
%   a,b      endpoints of an interval such that f(a)f(b) < 0
% Output       
%   x        vector of root approximations (last is best)

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

fa = f(a);  fb = f(b);
x = [];
w = b-a;                                % interval width
m = b;  fm = fb;                        % set up 1st iteration

while (w > xtol) & (abs(fm) > funtol)
  w = w/2;
  m = a + w;                            % midpoint of interval
  fm = f(m);
  % Update endpoints.
  if sign(fa) ~= sign(fm)
    b = m;
  else
    a = m;
    fa = fm;
  end
  x(end+1) = m;
end