[f,b,c,r,re,ae,iflag]=fzero(f,b,c,r,re,ae,iflag);

function [f,b,c,r,re,ae,iflag]=fzero(f,b,c,r,re,ae,iflag);
%***BEGIN PROLOGUE  FZERO
%***PURPOSE  Search for a zero of a function F(X) in a given interval
%            (B,C).  It is designed primarily for problems where F(B)
%            and F(C) have opposite signs.
%***LIBRARY   SLATEC
%***CATEGORY  F1B
%***TYPE      SINGLE PRECISION (FZERO-S, DFZERO-D)
%***KEYWORDS  BISECTION, NONLINEAR EQUATIONS, ROOTS, ZEROS
%***AUTHOR  Shampine, L. F., (SNLA)
%           Watts, H. A., (SNLA)
%***DESCRIPTION
%
%     FZERO searches for a zero of a REAL function F(X) between the
%     given REAL values B and C until the width of the interval (B,C)
%     has collapsed to within a tolerance specified by the stopping
%     criterion,
%        ABS(B-C) .LE. 2.*(RW*ABS(B)+AE).
%     The method used is an efficient combination of bisection and the
%     secant rule and is due to T. J. Dekker.
%
%     Description Of Arguments
%
%   F     :EXT   - Name of the REAL external function.  This name must
%                  be in an EXTERNAL statement in the calling program.
%                  F must be a function of one REAL argument.
%
%   B     :INOUT - One end of the REAL interval (B,C).  The value
%                  returned for B usually is the better approximation
%                  to a zero of F.
%
%   C     :INOUT - The other end of the REAL interval (B,C)
%
%   R     :OUT   - A (better) REAL guess of a zero of F which could help
%                  in speeding up convergence.  If F(B) and F(R) have
%                  opposite signs, a root will be found in the interval
%                  (B,R); if not, but F(R) and F(C) have opposite signs,
%                  a root will be found in the interval (R,C);
%                  otherwise, the interval (B,C) will be searched for a
%                  possible root.  When no better guess is known, it is
%                  recommended that r be set to B or C, since if R is
%                  not interior to the interval (B,C), it will be
%                  ignored.
%
%   RE    :IN    - Relative error used for RW in the stopping criterion.
%                  If the requested RE is less than machine precision,
%                  then RW is set to approximately machine precision.
%
%   AE    :IN    - Absolute error used in the stopping criterion.  If
%                  the given interval (B,C) contains the origin, then a
%                  nonzero value should be chosen for AE.
%
%   IFLAG :OUT   - A status code.  User must check IFLAG after each
%                  call.  Control returns to the user from FZERO in all
%                  cases.
%
%                1  B is within the requested tolerance of a zero.
%                   The interval (B,C) collapsed to the requested
%                   tolerance, the function changes sign in (B,C), and
%                   F(X) decreased in magnitude as (B,C) collapsed.
%
%                2  F(B) = 0.  However, the interval (B,C) may not have
%                   collapsed to the requested tolerance.
%
%                3  B may be near a singular point of F(X).
%                   The interval (B,C) collapsed to the requested tol-
%                   erance and the function changes sign in (B,C), but
%                   F(X) increased in magnitude as (B,C) collapsed, i.e.
%                     ABS(F(B out)) .GT. MAX(ABS(F(B in)),ABS(F(C in)))
%
%                4  No change in sign of F(X) was found although the
%                   interval (B,C) collapsed to the requested tolerance.
%                   The user must examine this case and decide whether
%                   B is near a local minimum of F(X), or B is near a
%                   zero of even multiplicity, or neither of these.
%
%                5  Too many (.GT. 500) function evaluations used.
%
%***REFERENCES  L. F. Shampine and H. A. Watts, FZERO, a root-solving
%                 code, Report SC-TM-70-631, Sandia Laboratories,
%                 September 1970.
%               T. J. Dekker, Finding a zero by means of successive
%                 linear interpolation, Constructive Aspects of the
%                 Fundamental Theorem of Algebra, edited by B. Dejon
%                 and P. Henrici, Wiley-Interscience, 1969.
%***ROUTINES CALLED  R1MACH
%***REVISION HISTORY  (YYMMDD)
%   700901  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  FZERO
persistent a acbs acmb aw cmb er fa fb fc fx fz gt ic kount p q rw t tol z ; 

if isempty(a), a=0; end;
if isempty(acbs), acbs=0; end;
if isempty(acmb), acmb=0; end;
if isempty(aw), aw=0; end;
if isempty(cmb), cmb=0; end;
if isempty(er), er=0; end;
if isempty(fa), fa=0; end;
if isempty(fb), fb=0; end;
if isempty(fc), fc=0; end;
if isempty(fx), fx=0; end;
if isempty(fz), fz=0; end;
if isempty(p), p=0; end;
if isempty(q), q=0; end;
if isempty(rw), rw=0; end;
if isempty(t), t=0; end;
if isempty(tol), tol=0; end;
if isempty(z), z=0; end;
if isempty(ic), ic=0; end;
if isempty(kount), kount=0; end;
if isempty(gt), gt=0; end;
%***FIRST EXECUTABLE STATEMENT  FZERO
%
%   ER is two times the computer unit roundoff value which is defined
%   here by the function R1MACH.
%
er = 2.0e0.*r1mach(4);
%
%   Initialize.
%
z = r;
if( r<=min(b,c) || r>=max(b,c) )
z = c;
end;
rw = max(re,er);
aw = max(ae,0.0e0);
ic = 0;
t = z;
fz = f(t);
fc = fz;
t = b;
fb = f(t);
kount = 2;
if( (abs(1.0e0).*sign(fz))~=(abs(1.0e0).*sign(fb)) )
c = z;
elseif( z~=c ) ;
t = c;
fc = f(t);
kount = 3;
if( (abs(1.0e0).*sign(fz))~=(abs(1.0e0).*sign(fc)) )
b = z;
fb = fz;
end;
end;
a = c;
fa = fc;
acbs = abs(b-c);
fx = max(abs(fb),abs(fc));
%
while( true );
if( abs(fc)<abs(fb) )
%
%   Perform interchange.
%
a = b;
fa = fb;
b = c;
fb = fc;
c = a;
fc = fa;
end;
%
cmb = 0.5e0.*(c-b);
acmb = abs(cmb);
tol = rw.*abs(b) + aw;
%
%   Test stopping criterion and function count.
%
if( acmb<=tol )
%
%   Finished.  Process results for proper setting of IFLAG.
%
if( (abs(1.0e0).*sign(fb))==(abs(1.0e0).*sign(fc)) )
iflag = 4;
elseif( abs(fb)>fx ) ;
iflag = 3;
else;
iflag = 1;
end;
return;
elseif( fb==0.0e0 ) ;
iflag = 2;
return;
else;
if( kount>=500 )
break;
end;
%
%   Calculate new iterate implicitly as B+P/Q, where we arrange
%   P .GE. 0.  The implicit form is used to prevent overflow.
%
p =(b-a).*fb;
q = fa - fb;
if( p<0.0e0 )
p = -p;
q = -q;
end;
%
%   Update A and check for satisfactory reduction in the size of the
%   bracketing interval.  If not, perform bisection.
%
a = b;
fa = fb;
ic = fix(ic + 1);
gt=0;
if( ic>=4 )
if( 8.0e0.*acmb>=acbs )
%
%   use bisection (C+B)/2.
%
b = b + cmb;
gt=1;
else;
ic = 0;
acbs = acmb;
end;
end;
%
%   Test for too small a change.
%
if(gt==0)
if( p<=abs(q).*tol )
%
%   Increment by TOLerance.
%
b = b + (abs(tol).*sign(cmb));
%
%   Root ought to be between B and (C+B)/2.
%
elseif( p>=cmb.*q ) ;
b = b + cmb;
else;
%
%   use secant rule.
%
b = b + p./q;
end;
end;
%
%   Have completed computation for new iterate B.
%
t = b;
fb = f(t);
kount = fix(kount + 1);
%
%   Decide whether next step is interpolation or extrapolation.
%
if( (abs(1.0e0).*sign(fb))==(abs(1.0e0).*sign(fc)) )
c = a;
fc = fa;
end;
end;
end;
iflag = 5;
return;
end %subroutine fzero
%DECK GAMIC