function [c w k] = bisect(a,b,f,M,epsilon,delta)
% BISECT implements bisection following pp 76-78 of Kincaid&Cheney
% usage examples:
%   >> bisect(0,1,@(x) cos(x)-x)
%   >> bisect(0,1,@(x) cos(x)-x,10,1.0e-6,1.0e-6)  % only allow 10 steps; ask for
%                                                  % 6 digit accuracy and residual
%   >> bisect(0,1,@fff)  % if function is in an m-file:
%                        %    function y = fff(x)
%                        %    y = cos(x) - x;

if a >= b, error('required: b should be greater than a'), end

u = feval(f,a);
v = feval(f,b);

if sign(u) == sign(v)
  error('required bracket: sign(f(a)) should not be same as sign(f(b))')
end

e = b - a;

if nargin < 4, M = 50; end
if nargin < 5, epsilon = 1.0e-14; end
if nargin < 6, delta = max((b-a) * 1.0e-14,1.0e-14); end

for k = 1:M
  e = e / 2.0;
  c = a + e;
  w = feval(f,c);
  if abs(e) < delta, return; end
  if abs(w) < epsilon, return; end
  if sign(w) ~= sign(u)
    b = c;
    v = w;
  else
    a = c;
    u = w;
  end
end

