function [c,yc,err,P] = bisect(f,a,b,delta)
% [c,yc,err] = bisect(f,a,b,delta)
% [c,yc,err,P] = bisect(f,a,b,delta)
% The bisection method is used to locate a root.
% f is the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% delta is the tolerance, input.
% c is the root, output.
% yc is the function value, output.
% err is the error estimate for c, output.
% P is the is the matrix of endpoint iterations, output.
P = [a b];
ya = feval(f,a);
yb = feval(f,b);
if ya*yb > 0, break, end
max1 = 1 + round((log(b-a)-log(delta))/log(2));
for k=1:max1,
c = (a+b)/2;
yc = feval(f,c);
if yc == 0,
a = c;
b = c;
elseif yb*yc > 0,
b = c;
yb = yc;
else
a = c;
ya = yc;
end
P = [P;a b];
if b-a < delta, break, end
end
c = (a+b)/2;
yc = feval(f,c);
err = abs(b-a)/2;