[x,y,yprime,neq,res,jac,h,wt,idid,rpar,ipar,phi,delta,e,wm,iwm,hmin,uround,nonneg,ntemp]=ddaini(x,y,yprime,neq,res,jac,h,wt,idid,rpar,ipar,phi,delta,e,wm,iwm,hmin,uround,nonneg,ntemp);

function [x,y,yprime,neq,res,jac,h,wt,idid,rpar,ipar,phi,delta,e,wm,iwm,hmin,uround,nonneg,ntemp]=ddaini(x,y,yprime,neq,res,jac,h,wt,idid,rpar,ipar,phi,delta,e,wm,iwm,hmin,uround,nonneg,ntemp);
%***BEGIN PROLOGUE  DDAINI
%***SUBSIDIARY
%***PURPOSE  Initialization routine for DDASSL.
%***LIBRARY   SLATEC (DASSL)
%***TYPE      doubleprecision (SDAINI-S, DDAINI-D)
%***AUTHOR  Petzold, Linda R., (LLNL)
%***DESCRIPTION
%-----------------------------------------------------------------
%     DDAINI TAKES ONE STEP OF SIZE H OR SMALLER
%     WITH THE BACKWARD EULER METHOD, TO
%     FIND YPRIME.  X AND Y ARE UPDATED TO BE CONSISTENT WITH THE
%     NEW STEP.  A MODIFIED DAMPED NEWTON ITERATION IS USED TO
%     SOLVE THE CORRECTOR ITERATION.
%
%     THE INITIAL GUESS FOR YPRIME IS USED IN THE
%     PREDICTION, AND IN FORMING THE ITERATION
%     MATRIX, BUT IS NOT INVOLVED IN THE
%     ERROR TEST. THIS MAY HAVE TROUBLE
%     CONVERGING IF THE INITIAL GUESS IS NO
%     GOOD, OR IF G(X,Y,YPRIME) DEPENDS
%     NONLINEARLY ON YPRIME.
%
%     THE PARAMETERS REPRESENT:
%     X --         INDEPENDENT VARIABLE
%     Y --         SOLUTION VECTOR AT X
%     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
%     NEQ --       NUMBER OF EQUATIONS
%     H --         STEPSIZE. IMDER MAY USE A STEPSIZE
%                  SMALLER THAN H.
%     WT --        VECTOR OF WEIGHTS FOR ERROR
%                  CRITERION
%     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS
%                  IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY
%                  IDID=-12 -- DDAINI FAILED TO FIND YPRIME
%     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS
%                  THAT ARE NOT ALTERED BY DDAINI
%     PHI --       WORK SPACE FOR DDAINI
%     DELTA,E --   WORK SPACE FOR DDAINI
%     WM,IWM --    REAL AND INTEGER ARRAYS STORING
%                  MATRIX INFORMATION
%
%-----------------------------------------------------------------
%***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV
%***REVISION HISTORY  (YYMMDD)
%   830315  DATE WRITTEN
%   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
%   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
%   901026  Added explicit declarations for all variables and minor
%           cosmetic changes to prologue.  (FNF)
%   901030  Minor corrections to declarations.  (FNF)
%***end PROLOGUE  DDAINI
%
persistent cj convgd damp delnrm err firstCall gt i ier ires jcalc lnje lnre m maxit mjac ncf nef nsf oldnrm r rate s xold ynorm ; if isempty(firstCall),firstCall=1;end; 

ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
iwm_shape=size(iwm);iwm=reshape(iwm,1,[]);
if isempty(gt), gt=0; end;
y_shape=size(y);y=reshape(y,1,[]);
yprime_shape=size(yprime);yprime=reshape(yprime,1,[]);
wt_shape=size(wt);wt=reshape(wt,1,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
phi_shape=size(phi);phi=reshape([phi(:).',zeros(1,ceil(numel(phi)./prod([neq])).*prod([neq])-numel(phi))],neq,[]);
delta_shape=size(delta);delta=reshape(delta,1,[]);
e_shape=size(e);e=reshape(e,1,[]);
wm_shape=size(wm);wm=reshape(wm,1,[]);
%
%
if isempty(i), i=0; end;
if isempty(ier), ier=0; end;
if isempty(ires), ires=0; end;
if isempty(jcalc), jcalc=0; end;
if isempty(m), m=0; end;
if isempty(maxit), maxit=0; end;
if isempty(mjac), mjac=0; end;
if isempty(ncf), ncf=0; end;
if isempty(nef), nef=0; end;
if isempty(nsf), nsf=0; end;
if isempty(cj), cj=0; end;
if isempty(damp), damp=0; end;
if isempty(delnrm), delnrm=0; end;
if isempty(err), err=0; end;
if isempty(oldnrm), oldnrm=0; end;
if isempty(r), r=0; end;
if isempty(rate), rate=0; end;
if isempty(s), s=0; end;
if isempty(xold), xold=0; end;
if isempty(ynorm), ynorm=0; end;
if isempty(convgd), convgd=false; end;
%
if isempty(lnre), lnre=12 ; end;
if isempty(lnje), lnje=13 ; end;
%
if firstCall,   maxit=[10];  end;
if firstCall,  mjac=[5];  end;
if firstCall,   damp=[0.75d0];  end;
firstCall=0;
%
%
%---------------------------------------------------
%     BLOCK 1.
%     INITIALIZATIONS.
%---------------------------------------------------
%
%***FIRST EXECUTABLE STATEMENT  DDAINI
idid = 1;
nef = 0;
ncf = 0;
nsf = 0;
xold = x;
[ynorm ,neq,y,wt,rpar,ipar]=ddanrm(neq,y,wt,rpar,ipar);
%
%     SAVE Y AND YPRIME IN PHI
for i = 1 : neq;
phi(i,1) = y(i);
phi(i,2) = yprime(i);
end; i = fix(neq+1);
%
%
%----------------------------------------------------
%     BLOCK 2.
%     DO ONE BACKWARD EULER STEP.
%----------------------------------------------------
%
%     SET UP FOR START OF CORRECTOR ITERATION
while( true );
cj = 1.0d0./h;
x = x + h;
%
%     PREDICT SOLUTION AND DERIVATIVE
for i = 1 : neq;
y(i) = y(i) + h.*yprime(i);
end; i = fix(neq+1);
%
jcalc = -1;
m = 0;
convgd = true;
%
%
%     CORRECTOR LOOP.
gt=0;
while( true );
iwm(lnre) = fix(iwm(lnre) + 1);
ires = 0;
%
[x,y,yprime,delta,ires,rpar,ipar]=res(x,y,yprime,delta,ires,rpar,ipar);
if( ires<0 )
%
%
%     EXITS FROM CORRECTOR LOOP.
convgd = false;
gt=1;
break;
else;
%
%
%     EVALUATE THE ITERATION MATRIX
if( jcalc==-1 )
iwm(lnje) = fix(iwm(lnje) + 1);
jcalc = 0;
[neq,x,y,yprime,delta,cj,h,ier,wt,e,wm,iwm,res,ires,uround,jac,rpar,ipar,ntemp]=ddajac(neq,x,y,yprime,delta,cj,h,ier,wt,e,wm,iwm,res,ires,uround,jac,rpar,ipar,ntemp);
%
s = 1000000.0d0;
if( ires<0 )
convgd = false;
gt=1;
break;
elseif( ier~=0 ) ;
convgd = false;
gt=1;
break;
else;
nsf = 0;
end;
end;
%
%
%
%     MULTIPLY RESIDUAL BY DAMPING FACTOR
for i = 1 : neq;
delta(i) = delta(i).*damp;
end; i = fix(neq+1);
%
%     COMPUTE A NEW ITERATE (BACK SUBSTITUTION)
%     STORE THE CORRECTION IN DELTA
%
[neq,delta,wm,iwm]=ddaslv(neq,delta,wm,iwm);
%
%     UPDATE Y AND YPRIME
for i = 1 : neq;
y(i) = y(i) - delta(i);
yprime(i) = yprime(i) - cj.*delta(i);
end; i = fix(neq+1);
%
%     TEST FOR CONVERGENCE OF THE ITERATION.
%
[delnrm ,neq,delta,wt,rpar,ipar]=ddanrm(neq,delta,wt,rpar,ipar);
if( delnrm<=100.0d0.*uround.*ynorm )
break;
end;
%
if( m>0 )
%
rate =(delnrm./oldnrm).^(1.0d0./m);
if( rate>0.90d0 )
convgd = false;
gt=1;
break;
else;
s = rate./(1.0d0-rate);
end;
else;
oldnrm = delnrm;
end;
%
if( s.*delnrm<=0.33d0 )
break;
end;
%
%
%     THE CORRECTOR HAS NOT YET CONVERGED. UPDATE
%     M AND AND TEST WHETHER THE MAXIMUM
%     NUMBER OF ITERATIONS HAVE BEEN TRIED.
%     EVERY MJAC ITERATIONS, GET A NEW
%     ITERATION MATRIX.
%
m = fix(m + 1);
if( m>=maxit )
convgd = false;
gt=1;
break;
else;
%
if(((fix(m./mjac)).*mjac)==m )
jcalc = -1;
end;
end;
end;
end;
%
%
%     THE ITERATION HAS CONVERGED.
%     CHECK NONNEGATIVITY CONSTRAINTS
if(gt==0)
if( nonneg~=0 )
for i = 1 : neq;
delta(i) = min(y(i),0.0d0);
end; i = fix(neq+1);
%
[delnrm ,neq,delta,wt,rpar,ipar]=ddanrm(neq,delta,wt,rpar,ipar);
if( delnrm>0.33d0 )
convgd = false;
else;
%
for i = 1 : neq;
y(i) = y(i) - delta(i);
yprime(i) = yprime(i) - cj.*delta(i);
end; i = fix(neq+1);
end;
end;
end;
if( convgd )
%
%
%
%-----------------------------------------------------
%     BLOCK 3.
%     THE CORRECTOR ITERATION CONVERGED.
%     DO ERROR TEST.
%-----------------------------------------------------
%
for i = 1 : neq;
e(i) = y(i) - phi(i,1);
end; i = fix(neq+1);
[err ,neq,e,wt,rpar,ipar]=ddanrm(neq,e,wt,rpar,ipar);
%
if( err<=1.0d0 )
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yprime_shape=zeros(yprime_shape);yprime_shape(:)=yprime(1:numel(yprime_shape));yprime=yprime_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
phi_shape=zeros(phi_shape);phi_shape(:)=phi(1:numel(phi_shape));phi=phi_shape;
delta_shape=zeros(delta_shape);delta_shape(:)=delta(1:numel(delta_shape));delta=delta_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
return;
end;
end;
%
%
%
%--------------------------------------------------------
%     BLOCK 4.
%     THE BACKWARD EULER STEP FAILED. RESTORE X, Y
%     AND YPRIME TO THEIR ORIGINAL VALUES.
%     REDUCE STEPSIZE AND TRY AGAIN, IF
%     POSSIBLE.
%---------------------------------------------------------
%
x = xold;
for i = 1 : neq;
y(i) = phi(i,1);
yprime(i) = phi(i,2);
end; i = fix(neq+1);
%
if( convgd )
%
nef = fix(nef + 1);
r = 0.90d0./(2.0d0.*err+0.0001d0);
r = max(0.1d0,min(0.5d0,r));
h = h.*r;
if( abs(h)<hmin || nef>=10 )
idid = -12;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yprime_shape=zeros(yprime_shape);yprime_shape(:)=yprime(1:numel(yprime_shape));yprime=yprime_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
phi_shape=zeros(phi_shape);phi_shape(:)=phi(1:numel(phi_shape));phi=phi_shape;
delta_shape=zeros(delta_shape);delta_shape(:)=delta(1:numel(delta_shape));delta=delta_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
return;
end;
elseif( ier~=0 ) ;
nsf = fix(nsf + 1);
h = h.*0.25d0;
if( nsf>=3 || abs(h)<hmin )
break;
end;
elseif( ires>-2 ) ;
ncf = fix(ncf + 1);
h = h.*0.25d0;
if( ncf>=10 || abs(h)<hmin )
idid = -12;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yprime_shape=zeros(yprime_shape);yprime_shape(:)=yprime(1:numel(yprime_shape));yprime=yprime_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
phi_shape=zeros(phi_shape);phi_shape(:)=phi(1:numel(phi_shape));phi=phi_shape;
delta_shape=zeros(delta_shape);delta_shape(:)=delta(1:numel(delta_shape));delta=delta_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
return;
end;
else;
idid = -12;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yprime_shape=zeros(yprime_shape);yprime_shape(:)=yprime(1:numel(yprime_shape));yprime=yprime_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
phi_shape=zeros(phi_shape);phi_shape(:)=phi(1:numel(phi_shape));phi=phi_shape;
delta_shape=zeros(delta_shape);delta_shape(:)=delta(1:numel(delta_shape));delta=delta_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
return;
end;
end;
idid = -12;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yprime_shape=zeros(yprime_shape);yprime_shape(:)=yprime(1:numel(yprime_shape));yprime=yprime_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
phi_shape=zeros(phi_shape);phi_shape(:)=phi(1:numel(phi_shape));phi=phi_shape;
delta_shape=zeros(delta_shape);delta_shape(:)=delta(1:numel(delta_shape));delta=delta_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
return;
%
%-------------end OF SUBROUTINE DDAINI----------------------
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yprime_shape=zeros(yprime_shape);yprime_shape(:)=yprime(1:numel(yprime_shape));yprime=yprime_shape;
wt_shape=zeros(wt_shape);wt_shape(:)=wt(1:numel(wt_shape));wt=wt_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
phi_shape=zeros(phi_shape);phi_shape(:)=phi(1:numel(phi_shape));phi=phi_shape;
delta_shape=zeros(delta_shape);delta_shape(:)=delta(1:numel(delta_shape));delta=delta_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
end %subroutine ddaini
%DECK DDAJAC