[ae,f,h,n,nq,iroot,re,t,yh,uround,b,c,fb,fc,y]=ddzro(ae,f,h,n,nq,iroot,re,t,yh,uround,b,c,fb,fc,y);

function [ae,f,h,n,nq,iroot,re,t,yh,uround,b,c,fb,fc,y]=ddzro(ae,f,h,n,nq,iroot,re,t,yh,uround,b,c,fb,fc,y);
%***BEGIN PROLOGUE  DDZRO
%***SUBSIDIARY
%***PURPOSE  DDZRO searches for a zero of a function F(N, T, Y, IROOT)
%            between the given 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).
%***LIBRARY   SLATEC (SDRIVE)
%***TYPE      doubleprecision (SDZRO-S, DDZRO-D, CDZRO-C)
%***AUTHOR  Kahaner, D. K., (NIST)
%             National Institute of Standards and Technology
%             Gaithersburg, MD  20899
%           Sutherland, C. D., (LANL)
%             Mail Stop D466
%             Los Alamos National Laboratory
%             Los Alamos, NM  87545
%***DESCRIPTION
%
%     This is a special purpose version of ZEROIN, modified for use with
%     the DDRIV package.
%
%     Sandia Mathematical Program Library
%     Mathematical Computing Services Division 5422
%     Sandia Laboratories
%     P. O. Box 5800
%     Albuquerque, New Mexico  87115
%     Control Data 6600 Version 4.5, 1 November 1971
%
%     PARAMETERS
%        F     - Name of the external function, which returns a
%                doubleprecision result.  This name must be in an
%                EXTERNAL statement in the calling program.
%        B     - One end of the interval (B, C).  The value returned for
%                B usually is the better approximation to a zero of F.
%        C     - The other end of the interval (B, C).
%        RE    - 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    - 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.
%
%***REFERENCES  L. F. Shampine and H. A. Watts, ZEROIN, a root-solving
%                 routine, SC-TM-70-631, Sept 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, 1969.
%***ROUTINES CALLED  DDNTP
%***REVISION HISTORY  (YYMMDD)
%   790601  DATE WRITTEN
%   900329  Initial submission to SLATEC.
%***end PROLOGUE  DDZRO
persistent a acbs acmb cmb er fa gt ic kount p q rw tol ; 

if isempty(ic), ic=0; end;
if isempty(kount), kount=0; end;
if isempty(gt), gt=0; end;
if isempty(a), a=0; end;
if isempty(acbs), acbs=0; end;
if isempty(acmb), acmb=0; end;
if isempty(cmb), cmb=0; end;
if isempty(er), er=0; end;
if isempty(fa), fa=0; end;
if isempty(p), p=0; end;
if isempty(q), q=0; end;
if isempty(rw), rw=0; end;
if isempty(tol), tol=0; end;
y_shape=size(y);y=reshape(y,1,[]);
yh_shape=size(yh);yh=reshape([yh(:).',zeros(1,ceil(numel(yh)./prod([n])).*prod([n])-numel(yh))],n,[]);
%***FIRST EXECUTABLE STATEMENT  DDZRO
er = 4.0d0.*uround;
rw = max(re,er);
ic = 0;
acbs = abs(b-c);
a = c;
fa = fc;
kount = 0;
%                                                    Perform interchange
while( true );
if( abs(fc)<abs(fb) )
a = b;
fa = fb;
b = c;
fb = fc;
c = a;
fc = fa;
end;
cmb = 0.5d0.*(c-b);
acmb = abs(cmb);
tol = rw.*abs(b) + ae;
%                                                Test stopping criterion
if( acmb<=tol )
break;
end;
if( kount>50 )
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.0d0 )
p = -p;
q = -q;
end;
%                          Update A and check for satisfactory reduction
%                          in the size of our bounding interval.
a = b;
fa = fb;
ic = fix(ic + 1);
gt=0;
if( ic>=4 )
if( 8.0d0.*acmb>=acbs )
%                                                                 Bisect
b = 0.5d0.*(c+b);
gt=1;
else;
ic = 0;
end;
end;
if(gt==0)
acbs = acmb;
%                                            Test for too small a change
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 ) ;
%                                                            Interpolate
b = b + p./q;
else;
%                                                                 Bisect
b = 0.5d0.*(c+b);
end;
end;
%                                             Have completed computation
%                                             for new iterate B.
[h,dumvar2,n,nq,t,b,yh,y]=ddntp(h,0,n,nq,t,b,yh,y);
fb = f(n,b,y,iroot);
if( n==0 )
break;
end;
if( fb==0.0d0 )
break;
end;
kount = fix(kount + 1);
%
%             Decide whether next step is interpolation or extrapolation
%
if( (abs(1.0d0).*sign(fb))==(abs(1.0d0).*sign(fc)) )
c = a;
fc = fa;
end;
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yh_shape=zeros(yh_shape);yh_shape(:)=yh(1:numel(yh_shape));yh=yh_shape;
end %subroutine ddzro
%DECK DE1