38 lines
1.0 KiB
Fortran
38 lines
1.0 KiB
Fortran
SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
|
|
INTEGER n,NMAX
|
|
REAL eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n)
|
|
EXTERNAL derivs
|
|
PARAMETER (NMAX=50)
|
|
CU USES derivs,rkck
|
|
INTEGER i
|
|
REAL errmax,h,htemp,xnew,yerr(NMAX),ytemp(NMAX),SAFETY,PGROW,
|
|
*PSHRNK,ERRCON
|
|
PARAMETER (SAFETY=0.9,PGROW=-.2,PSHRNK=-.25,ERRCON=1.89e-4)
|
|
h=htry
|
|
1 call rkck(y,dydx,n,x,h,ytemp,yerr,derivs)
|
|
errmax=0.
|
|
do 11 i=1,n
|
|
errmax=max(errmax,abs(yerr(i)/yscal(i)))
|
|
11 continue
|
|
errmax=errmax/eps
|
|
if(errmax.gt.1.)then
|
|
htemp=SAFETY*h*(errmax**PSHRNK)
|
|
h=sign(max(abs(htemp),0.1*abs(h)),h)
|
|
xnew=x+h
|
|
if(xnew.eq.x)pause 'stepsize underflow in rkqs'
|
|
goto 1
|
|
else
|
|
if(errmax.gt.ERRCON)then
|
|
hnext=SAFETY*h*(errmax**PGROW)
|
|
else
|
|
hnext=5.*h
|
|
endif
|
|
hdid=h
|
|
x=x+h
|
|
do 12 i=1,n
|
|
y(i)=ytemp(i)
|
|
12 continue
|
|
return
|
|
endif
|
|
END
|