SUBROUTINE ODEINT(YSTART,NVAR,X1,X2,EPS,H1,HMIN,NOK,NBAD,
     1          DERIVS,RKQS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c
c Integrator driver with adaptive stepsize control.  Integrate the NVAR 
c starting values YSTART from X1 to X2 with accuracy EPS, storing 
c intermediate results in the common block /PATH/.  H1 should be set as a 
c guessed first step size, HMIN as the minimum allowed step size (zero okay).
c On output, NOK and NBAD are the number of good and bad (but retried and 
c fixed) steps taken, and YSTART is replaced by values at the end of the
c integration interval.  DERIVS is the user-supplied subroutine for 
c calculating the right hand side derivative, while RKQS is the name of
c the stepper routine to be used.  PATH contains its own information about
c how often an intermediate value is to be stored.
c
      INTEGER nbad,nok,nvar,KMAXX,MAXSTP,NMAX
      DOUBLE PRECISION eps,h1,hmin,x1,x2,ystart(NVAR),TINY
      EXTERNAL DERIVS,RKQS
      PARAMETER (MAXSTP=10000,NMAX=50,KMAXX=20000,TINY=1.d-50)
      INTEGER i,kmax,kount,nstp
      DOUBLE PRECISION dxsav,h,hdid,hnext,x,xsav,dydx(NMAX),xp(KMAXX),
     1  y(NMAX),yp(NMAX,KMAXX),yscal(NMAX)
c user storage for intermediate results.
c remember to preset DXSAV and KMAX
      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.d0*dxsav
      DO 16 NSTP=1,MAXSTP
        CALL DERIVS(X,Y,DYDX)
C scaling used to monitor accuracy.  This general purpose choice can
c be modified if need be.
        DO 12 I=1,NVAR
          YSCAL(I)=ABS(Y(I))+ABS(H*DYDX(I))+TINY
12      CONTINUE
C store intermediate results
        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
c if step can overshoot end, cut down step size
        IF((X+H-X2)*(X+H-X1).GT.0.d0) 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
C are we done?
        IF ((X-X2)*(X2-X1).GE.0.d0) THEN
          DO 14 I=1,NVAR
            YSTART(I)=Y(I)
14        CONTINUE
          IF (KMAX.NE.0) THEN
C  save final step
            KOUNT=KOUNT+1
            XP(KOUNT)=X
            DO 15 I=1,NVAR
              YP(I,KOUNT)=Y(I)
15          CONTINUE
          ENDIF
C    normal return
          RETURN
        ENDIF
        IF (ABS(HNEXT).LT.HMIN) HNEXT=HMIN
        H=HNEXT
16    CONTINUE
      PAUSE 'Too many steps!'
      RETURN
      END
c