      SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,
     *rkqs)
      INTEGER nbad,nok,nvar,KMAXX,MAXSTP,NMAX
      REAL eps,h1,hmin,x1,x2,ystart(nvar),TINY
      EXTERNAL derivs,rkqs
      PARAMETER (MAXSTP=10000,NMAX=50,KMAXX=200,TINY=1.e-30)
C	Runge-Kutta driver with adaptive stepsize control.  Integrate the 
C	starting values ystart(1:nvar)from x1 to x2 with accuracy eps, 
C	storing intermediate results in the common block /path/
C	h1 should be set as a guessed first stepsize.  Hmin as the minimum 
C	allowed stesize (can be zero)  On output nok and nbad are the number 
C	of good and bad (but retried and fixed) steps taken, and ystart is 
C	replaced by values at the end of the integration interval.  Derivs 
C	is the user-supplied subrouting for calculating the right hand side 
C	derivative, while RKQS is the name of the stepper routine to be 
C	used. /path/ contains its own information about how often an 
C	intermediate value is to be stored.  
      INTEGER i,kmax,kount,nstp
      REAL dxsav,h,hdid,hnext,x,xsav,dydx(NMAX),xp(KMAXX),y(NMAX),
     *yp(NMAX,KMAXX),yscal(NMAX)
      COMMON /path/ kmax,kount,dxsav,xp,yp
      x=x1
      h=sign(h1,x2-x1)
      nok=0
      nbad=0
      kount=0
      do 11 i=1,nvar
        y(i)=ystart(i)
11    continue
      if (kmax.gt.0) xsav=x-2.*dxsav
      do 16 nstp=1,MAXSTP
        call derivs(x,y,dydx)
        do 12 i=1,nvar
          yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY
12      continue
        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 13 i=1,nvar
                yp(i,kount)=y(i)
13            continue
              xsav=x
            endif
          endif
        endif
        if((x+h-x2)*(x+h-x1).gt.0.) h=x2-x
        call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
        if(hdid.eq.h)then
          nok=nok+1
        else
          nbad=nbad+1
        endif
        if((x-x2)*(x2-x1).ge.0.)then
          do 14 i=1,nvar
            ystart(i)=y(i)
14        continue
          if(kmax.ne.0)then
            kount=kount+1
            xp(kount)=x
            do 15 i=1,nvar
              yp(i,kount)=y(i)
15          continue
          endif
          return
        endif
        if(abs(hnext).lt.hmin) THEN
	  PRINT*,'PAUSE AT X=',X,' HNEXT, HMIN',HNEXT,HMIN
	  pause 'stepsize smaller than minimum in odeint'
	ENDIF
        h=hnext
16    continue
      pause 'too many steps in odeint'
      return
      END
