function x_star = bisect( func, x_lft, x_rgt, ftol )
%
%   Find a solution to the equation f(x) = 0 for the function
%   with handle func using the method of bisection.  
%   Begin with the interval [ x_lft, x_rgt ] and return the solution 
%   with relative error ftol (optional; default 1e-3 takes 10 iterations).

if !exist('ftol','var')
   ftol = 1e-3;  % default fractional tolerance, ftol
end

f_lft = func( x_lft );    % f_lft = f(x_lft)
f_rgt = func( x_rgt );    % f_rgt = f(x_rgt)

if f_lft*f_rgt > 0          % ooops.
   return                   % bail.
end

while abs( x_rgt - x_lft ) > ftol*max( abs( [ x_lft x_rgt ] ) )
   x_mp = 0.5*( x_lft + x_rgt );
   f_mp = func( x_mp );
   
   if f_mp == 0           % got lucky!
      x_star = x_mp;
      return
   end
   
   if f_mp*f_rgt > 0
      x_rgt = x_mp;
      f_rgt = f_mp;
   else
      x_lft= x_mp;
      f_lft= f_mp;
   end
end

x_star = 0.5*( x_lft + x_rgt );
