      program vdpol
      INTEGER NSTEP,NVAR
      INTEGER KMAXX,NMAX
      PARAMETER(NVAR=2,NSTEP=6000,KMAXX=9000,NMAX=50)
      INTEGER i,j,kmax,kount,nbad,nok,nrhs
      REAL dxsav,eps,h1,hmin,x1,x2,x,y
c     REAL x(NSTEP),y(NVAR,NSTEP),
      REAL vstart(NVAR),time,dt,tmp
      COMMON /path/ kmax,kount,dxsav,x(KMAXX),y(NMAX,KMAXX)
      COMMON nrhs
      EXTERNAL derivs, rkqs
 
      nrhs=0
      x1=0.0
      x2=3000.0
      vstart(1)=2.0
      vstart(2)=0.0

      eps=1.0e-6
      h1=.1
      hmin=0.0
      kmax=KMAXX
      dxsav=(x2-x1)/(2.0*x2)

c      call rkdumb(vstart,NVAR,x1,x2,NSTEP,derivs)
      call odeint(vstart,NVAR,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs)
      write(*,'(/1x,a,t30,i9)') 'Successful steps:',nok
      write(*,'(1x,a,t30,i9)') 'Bad steps:',nbad
     
      write(*,'(1x,a,t30,i6)') 'Stored intermediate values:',kount
      i=nrhs
      write(*,'(1x,a,t30,i9)') 'Function evaluations:', i

      open (unit=10,file="vdpolx.dat",status="unknown")
      open (unit=11,file="vdpoly.dat",status="unknown")
      time=0
      tmp=nstep
c      dt=(x2-x1)/tmp
c      write(10,100) time,vstart(1)
c      write(11,100) time,vstart(2)
      do 1 i=1,kount
c         time=time+dt
      write(10,100) x(i),y(1,i)
      write(11,100) x(i),y(2,i)
1     continue

      close(10)
      close(11)
 100  Format (2(2x,f10.5))   
      stop
      end

      subroutine derivs(x,y,dydx)
      INTEGER nrhs
      real x,y(*),dydx(*)
      COMMON nrhs
      nrhs=nrhs+1
      dydx(1)=y(2)
      dydx(2)=1000*(1-y(1)**2)*y(2)-y(1)
      return
      end


      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=1500000,NMAX=50,KMAXX=9000,TINY=1.e-30)
      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) pause
     *'stepsize smaller than minimum in odeint'
        h=hnext
16    continue
      pause 'too many steps in odeint'
      return
      END
