function x = bisec(f, a, b)
% BISEC   bisection method for a scalar equation f(x) = 0
%
% Input:
%   f        objective function 
%   a, b     initial root approximations, a < b
%
% Output       
%   x        vector of root approximations (last is best)

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

    x = [a b];
    dx = b - a;
    fa = f(a);
    fb = f(b);
    fevals = 2;
    fc = Inf;

    while (dx > xtol) && (abs(fc) > funtol)
      c = (a + b)/2;
      fc = f(c);
      fevals = fevals + 1;
      x(fevals) = c;
      if fc == 0
          break
      end
      if sign(fc)*sign(fa) < 0
          b = c;
          fb = fc;
      else
          a = c;
          fa = fc;
      end
      dx = b - a;
    end
end