      program vdpol
      INTEGER NSTEP,NVAR
      INTEGER KMAXX,NMAX
      PARAMETER(NVAR=2,NSTEP=6000,KMAXX=200,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=20.0
      vstart(1)=2.0
      vstart(2)=0.0

      eps=1.0e-5
      h1=.1
      hmin=0.0
      kmax=200
      dxsav=(x2-x1)/NSTEP

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,i3)') 'Successful steps:',nok
      write(*,'(1x,a,t30,i3)') 'Bad steps:',nbad
      write(*,'(1x,a,t30,i3)') 'Function evaluations:',nrhs
      write(*,'(1x,a,t30,i3)') 'Stored intermediate values:',kount

      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)

      real x,y(*),dydx(*)

      dydx(1)=y(2)
      dydx(2)=(1-y(1)**2)*y(2)-y(1)
    

      return
      end


