subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs,ierror)
INTEGER nbad,nok,nvar
double precision eps,h1,hmin,x1,x2,ystart(nvar)
INTEGER i,nstp
double precision h,hdid,hnext,x,xsav,dydx(NMAX),y(NMAX),yscal(NMAX)
integer kmax, kount
double precision dxsav, yp(NMAX,KMAXX), xp(KMAXX)
! Can be 0: success, 1: hmin too small, 2: MAXSTP exceeded
integer, intent(out) :: ierror
EXTERNAL derivs,rkqs
COMMON /path/ kmax,kount,dxsav,xp,yp
x = x1
h = sign(h1,x2-x1)
nok = 0
nbad = 0
kount = 0
xsav = 0.d0
do i=1,nvar
y(i)=ystart(i)
end do
if (kmax.gt.0) xsav=x-2.d0*dxsav
do nstp=1,MAXSTP
call derivs(x,y,dydx)
do i=1,nvar
yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY
end do
if(kmax.gt.0)then
if(abs(x-xsav).gt.abs(dxsav)) then
if(kount.lt.kmax-1)then
kount=kount+1
xp(kount)=x
do i=1,nvar
yp(i,kount)=y(i)
end do
xsav=x
end if
end if
end if
if((x+h-x2)*(x+h-x1).gt.0.d0) h=x2-x
if (any(y(1:nvar) .ne. y(1:nvar)) .or. any(dydx(1:nvar) .ne. &
dydx(1:nvar)) .or. x .ne. x) then
write(*,*) "ODEINT BEFORE CALL RKQS: NaN in x, dydx, or y!"
write(*,*) "x",x
write(*,*) "y",y
write(*,*) "dydx",dydx
write(*,*) "exiting..."
ierror = 2
return
end if
call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
if (any(y(1:nvar) .ne. y(1:nvar)) .or. any(dydx(1:nvar) .ne. &
dydx(1:nvar)) .or. x .ne. x) then
write(*,*) "ODEINT AFTER CALL RKQS: NaN in x, dydx, or y!"
write(*,*) "x",x
write(*,*) "y",y
write(*,*) "dydx",dydx
write(*,*) "exiting..."
ierror = 2
return
end if
if(hdid.eq.h)then
nok=nok+1
else
nbad=nbad+1
end if
if((x-x2)*(x2-x1).ge.0.d0)then
do i=1,nvar
ystart(i)=y(i)
end do
if(kmax.ne.0)then
kount=kount+1
xp(kount)=x
do i=1,nvar
yp(i,kount)=y(i)
end do
end if
ierror = 0
return
end if
if(abs(hnext).lt.hmin) then
ierror = 1
end if
h=hnext
end do
ierror = 2
end subroutine odeint