%
% Find roots if f(x) = 0 using the Bisection Method
%
function [c,yc,err,H] = bisect(f,a,b,delta)

% Inputs
%   f       name of the function
%   a       left endpoint
%   b       right endpoint
%   delta   convergence tolerance
% Output
%   c       solution: the root
%   yc      solution: the function value
%   err     error estimate in c
%   H       History vector of the iterations


ya = feval(f,a); yb = feval(f,b);
H = [a b ya];
if ya*yb > 0, end
maxit = 1 + round((log(b-a)-log(delta))/log(2));
for it = 1:maxit
  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
  H = [H;a b yc];
  if b-a < delta, break, end
end
c  = (a+b)/2; yc = feval(f,c);
err = abs(b-a)/2;

