function [b,c,residual,flag] = Zero(f,b,c,abserr,relerr)
%  Zero computes a root of the nonlinear equation f(x) = 0 when f(x)
%  is a continuous real function of a single real variable x.  The
%  method used is a combination of bisection and the secant rule.
%
%    Input parameters:
%       f     = the name of the function f(x).
%       b,c   = values of x such that f(b)*f(c) <= 0.
%       abserr,relerr = absolute and relative error tolerances.
%               The stopping criterion is:
%                 abs(b-c) <= 2*max(abserr,abs(b)*relerr).
%    Output parameters:
%       b,c       = see flag options.
%       residual  = value of final residual f(b).
%       flag      = 0 for normal return; f(b)*f(c) < 0 and the
%                     stopping criterion is met (or f(b) = 0).  b is
%                     always selected so that abs(f(b)) <= abs(f(c)).
%                 = 1 if too many function evaluations were made; in
%                     this version 500 are allowed.
%                 = 2 if abs(f(b)) is more than 100 times the larger
%                     of the residuals for the input b and c.  b and c
%                     probably approximate a pole rather than a root.


maxfe = 500;
false = 0;
true = ~false;
    
%  Test the input tolerances.
if relerr < 10*eps
  error('relerr < 10*eps is not allowed in Zero.')
end
if abserr < 0
  error('abserr < 0 is not allowed in Zero.')
end

%  Initialization.
count = 0;
width = abs(b - c);
a = c;  
fa = feval(f,a);
nfe = 1;
if fa == 0
  b = a;
  residual = 0;
  flag = 0;
  return;
end
fb = feval(f,b);
nfe = 2;
if fb == 0
  residual = 0;
  flag = 0;
  return;
end
if (fa > 0) == (fb > 0)
  error('Initial interval [b,c] does not bracket a root.')
end
initial_residual = max(abs(fa),abs(fb));
fc = fa;
while true

  if abs(fc) < abs(fb)

    %  Interchange b and c so that abs(f(b)) <= abs(f(c)).
    a = b;  fa = fb;
    b = c;  fb = fc;
    c = a;  fc = fa;
  end
  cmb = 0.5*(c-b);
  acmb = abs(cmb);
  tol = max(abserr, abs(b)*relerr);

  %  Test the stopping criterion and function count.
  if acmb <= tol
    residual = fb;
    if abs(residual) > 100*initial_residual
      flag = 2;
      return;
    else
      flag = 0;
      return;
    end
  end
  if nfe >= maxfe
    flag = 1;
    return;
  end

  %  Calculate new iterate implicitly as b+p/q where we arrange
  %  p >= 0.  The implicit form is used to prevent overflow.
  p = (b-a)*fb;
  q = fa-fb;
  if p < 0
    p = -p;  q = -q;
  end

  %  Update a; check if reduction in the size of bracketing
  %  interval is being reduced at a reasonable rate.  If not,
  %  bisect until it is.
  a = b;  fa = fb;
  count = count + 1;
  bisect = false;
  if count >= 4
    if 8*acmb >= width
      bisect = true;
    else
      count = 0;
      width = acmb;
    end
  end

  %  Test for too small a change.
  if p <= abs(q)*tol
    %  If the change is too small, increment by tol.
    b = b + tol * sign(cmb);
  else
    %  Root ought to be between b and (c+b)/2.
    if p < cmb*q
      b = b + p/q;      %    Use secant rule.   
    else
      bisect = true;    %    Use bisection.     
    end
  end

  if bisect
    b = c - cmb;
  end

  %  The computation for new iterate b has been completed.
  fb = feval(f,b);
  nfe = nfe + 1;
  if fb == 0
    c = b;
    residual = fb;
    flag = 0;
    return;
  end
  if (fb > 0) == (fc > 0)
    c = a;  fc = fa;
  end
end