subroutine fcn(n, x, y, yprime) 
      integer n
      double precision x, y(n), yprime(n)

      alpha     = 0.01
      drabbitdt =  +2
      dfoxdt    =  -1

      nrabbits = y(1)
      nfoxes   = y(2)

      yprime(1) = drabbitdt*nrabbits - alpha*nrabbits*nfoxes
      yprime(2) = dfoxdt*nfoxes  + alpha*nrabbits*nfoxes

      return
      end
      

      program preditorprey

      implicit none
      external fcn
      double precision time, timeend, endstep, dtime 
      integer step, nsteps
      integer neqns
      integer status
      parameter (neqns = 2)
      double precision population(neqns)
      double precision tol

c
c     this is just workspace for the integrator; 
c     fortran didn't have allocation of memory in the olden days.
c      
      integer wsize
      parameter (wsize = neqns)
      double precision workarray(wsize,9)
      double precision comm(24)

c
c     set initial conditions
c
      population(1) = 300
      population(2) = 150
       
c
c     set tolerance
c
      tol = 1.d-6
c
c     how long do we want to run this thing, and output how many
c     points?
c
      time = 0.
      timeend = 40.
      nsteps  = 200
      dtime = (timeend - time)/(1.d0*nsteps)

      print *, '#nstep  time  rabbits  foxes'
      print *, 0, time, population(1), population(2)
      status = 1
      do step=1, nsteps
        endstep = time + dtime
        call dverk (neqns, fcn, time, population, endstep, tol, 
     1               status, comm, wsize, workarray)
        if (status .ne. 3) then
           if (status .gt. 0) then 
              print *, '# interrupt occured -- wierd but probably ok '
           else  
              print *, '# status is ', status, '; crash!!!'
              stop
           endif
        endif
        print *, step, endstep, population(1), population(2)
      enddo 
      end