function b = froot(F,a,b,varargin)
%FROOT Simplified version of FZERO.
%   X=FROOT(F,A,B) tries to find a zero of F(x) between A and B such that F(A) 
%   and F(B) must have opposite signs. More precisely, FROOT returns one
%   endpoint of a small subinterval of [A,B] where F changes sign.
%
%   X=FROOT(F,A,B,P1,P2,...) passes parameters P1,P2,... to F.
%
%   Adapted version of FZEROTX from "Numerical Computing with MATLAB" by 
%   C. Moler.

fa = F(a,varargin{:});
fb = F(b,varargin{:});
if sign(fa) == sign(fb)
   error('Function must change sign on the interval.')
end
c = a;
fc = fa;
d = b - c;
e = d;

while fb ~= 0
   
   % The three current points a, b, and c, satisfy:
   %
   % - f(x) changes sign between a and b;
   % - abs(f(b)) <= abs(f(a)); and
   % - c equals previous b, so c might be equal to a.
   %
   % The next point is chosen from:
   %
   % - bisection point, (a+b)/2; or
   % - secant point determined by b and c; or
   % - inverse quadratic interpolation point determined by a, b, and c 
   %   if they are distinct.
   
   if sign(fa) == sign(fb)
      a = c; fa = fc;
      d = b - c; e = d;
   end

   if abs(fa) < abs(fb)
      c = b; b = a; a = c;
      fc = fb; fb = fa; fa = fc;
   end
   
   % Convergence test.
   m = .5*(a - b);
   tol = 2*eps*max(abs(b),1);
   if (abs(m) <= tol) || (fb == 0.0)
      break
   end
   
   % Choose bisection or interpolation.
   if (abs(e) < tol) || (abs(fc) <= abs(fb))  % Bisection.
      d = m;
      e = m;
   else                                       % Interpolation.
      s = fb/fc;
      if a == c % Linear interpolation (secant).
         p = 2*m*s;
         q = 1 - s;
      else      % Inverse quadratic interpolation.
         q = fc/fa;
         r = fb/fa;
         p = s*(2*m*q*(q - r) - (b - c)*(r - 1));
         q = (q - 1)*(r - 1)*(s - 1);
      end;
      if p > 0
         q = -q; 
      else
         p = -p; 
      end;
      % Is interpolated point acceptable?
      if (2*p < 3*m*q - abs(tol*q)) && (p < abs(.5*e*q))
         e = d;
         d = p/q;
      else
         d = m;
         e = m;
      end;
   end
   % Next point.
   c = b;
   fc = fb;
   if abs(d) > tol
      b = b + d;
   else
      b = b - sign(b-a)*tol;
   end
   fb = F(b,varargin{:});
   
en