C
C  This is the FORTRAN driver to generate a discrete solution and
C  a piecewise polynomial solution from DVERK for the predator-prey
C  problem.
C
C  The output from this program is a file that is a Matlab script
C  for generating the solution vector on a fine mesh suitable
C  for plotting in Matlab.
C
      integer i, j, k, ind, n, ref, nw
      double precision TT(200), SD(2,200), y(2), t, tend, c(30), tol,
     +    W(2,30), TTF(2000), SDF(2, 2000), delta
      external fcn
      data c/30*0.d0/
      tol = 1.d-2
      n = 2
      nw = 2
      ref = 20
      y(1) = 30.d0
      y(2) = 20.d0
      t = 0.d0
      tend = 40.d0
      TT(1) = t
      TTF(1) = t
      SD(1,1) = y(1)
      SD(2,1) = y(2)
      SDF(1,1) = SD(1,1)
      SDF(2,1) = SD(2, 1)
      print 110, 1, TTF(1), 1, SDF(1,1), SDF(2,1)
      print 130, 1, TT(1), 1, SD(1,1), SD(2,1)
      ind = 2
      c(9) = 1
      c(6) = 10.d0
      i = 1
      j = 1
 100  call dverk(n, fcn, t, y, tend, tol, ind, c, nw, W)
      if ( ind .eq. 6 ) go to 100
      if ((ind .eq. 3) .or. (ind .lt. 0)) go to 400
C
C The step is to be accepted (ind=5) so we store the discrete solution
C and compute (using intrp) the off-mesh interpolated solution on a fine
C grid of 'ref' equally spaced points.
C
         TT(i+1) = c(17)
         SD(1, i+1) = W(1,13)
         SD(2, i+1) = W(2,13)
         delta = c(18)/dfloat(ref)
         do 120 k = 1,ref
            j = (i-1)*ref + k + 1
            TTF(j) = t + delta*k
            call intrp(n, t, y, TTF(j), SDF(1,j), c(18), nw, W)
            print 110, j, TTF(j),j, SDF(1,j), SDF(2,j)
 110        format(1x, 'TTF(', i4, ') = ', E15.8, ';'/,
     +        1x, 'SDF( 1:2 ,', i4, ') = [', 2E15.8, "]';")
 120     continue
         i = i+1
         print 130, i, TT(i), i, SD(1,i), SD(2,i)
 130     format(1x, 'TT(', i3, ') = ', E15.8, ';'/,
     +     1x, 'SD( 1:2,', i3, ') = [', 2E15.8, "]';" )
         go to 100
 400  print 410, i, j
 410  format(1x, 'II = ', i4, ';' /1x, 'JJ = ', i4, ';' )
      stop
      end

      subroutine fcn(n, t, y, yp)
      double precision t, y(2), yp(2)
      yp(1) = y(1) - .1d0*y(1)*y(2) + .02d0*t
      yp(2) = -y(2) + .02d0*y(1)*y(2) + .008d0*t
      return
      end