      program vdpol
      INTEGER NSTEP,NVAR
      PARAMETER(NVAR=2,NSTEP=6000)
      INTEGER i,j
      REAL xx,x1,x2,y,vstart(NVAR),time,dt,tmp
      COMMON /path/ xx(NSTEP+2),y(NVAR,NSTEP+2)
      EXTERNAL derivs
 
      x1=0.0
      x2=20.0
      vstart(1)=2.0
      vstart(2)=0.0


      call rkdumb(vstart,NVAR,x1,x2,NSTEP,derivs)

      open (unit=10,file="vdpolx.dat",status="unknown")
      open (unit=11,file="vdpoly.dat",status="unknown")
      time=0
      tmp=nstep
      dt=(x2-x1)/tmp
c      write(10,100) time,vstart(1)
c      write(11,100) time,vstart(2)
      do 1 i=1,nstep
c         time=time+dt
      write(10,100) xx(i),y(1,i)
      write(11,100) xx(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


      SUBROUTINE rkdumb(vstart,nvar,x1,x2,nstep,derivs)
      INTEGER nstep,nvar,NMAX,NSTPMX
      PARAMETER (NMAX=2,NSTPMX=6000)
      REAL x1,x2,vstart(nvar),xx,y
      EXTERNAL derivs
      COMMON /path/ xx(NSTPMX+2),y(NMAX,NSTPMX+2)
CU    USES rk4
      INTEGER i,k
      REAL h,x,dv(NMAX),v(NMAX)
      do 11 i=1,nvar
        v(i)=vstart(i)
        y(i,1)=v(i)
11    continue
      xx(1)=x1
      x=x1
      h=(x2-x1)/nstep
      do 13 k=1,nstep
        call derivs(x,v,dv)
        call rk4(v,dv,nvar,x,h,v,derivs)
        if(x+h.eq.x)pause 'stepsize not significant in rkdumb'
        x=x+h
        xx(k+1)=x
        do 12 i=1,nvar
          y(i,k+1)=v(i)
12      continue
13    continue
      return
      END
