c   example program for epsode -- robertson 3-species chemical
c     kinetics problem
c     taken from
c        hindmarsh and byrne, applications of episode: an experimental
c           package for the integration of systems of ordinary differ-
c           ential equations, from numerical methods for differential
c           systems, academic press 1976.
c
      implicit real*8 (a-h,o-z)
      dimension y0(3)
      common/epcom9/ hused, nqused, nstep, nfe, nje
      external diffun, jacob
c
      n = 3
      eps = 2.0d-7
      ierror = 1
      do 70 miter = 1, 2
         mf = 20 + miter
         t0 = 0.0d0
         tout = 0.4d0
         h0 = 1.0d-9
         y0(1) = 1.0d0
         y0(2) = 0.0d0
         y0(3) = 0.0d0
         write(6,1000) mf, eps
         nout = 12
         index = 1
         do 60 iout = 1, nout
            call epsode (diffun, jacob, n, t0, h0, y0, tout, eps,
     *                   ierror, mf, index)
            sy = y0(1) + y0(2) + y0(3) - 1.0d0
            write(6,1001) t0, nstep, nfe, nje, nqused, h0, y0(1),
     *                    y0(2), y0(3), sy
            if (dabs(sy) .gt. eps) stop
            if (index .ne. 0) stop
            tout = 10.0d0 * tout
   60    continue
   70 continue
      stop
 1000 format('1  robertson 3-species chemical kinetics problem'/
     *       '0mf =',i3,'   eps =',1pd10.1/
     *       '-   t',6x,'nstep',3x,'nfe',3x,'nje',2x,'nq',4x,'h',10x,
     *       'y1',10x,'y2',10x,'y3',5x,' sum(y)-1'//)
 1001 format(1x,1pd9.1,3i6,i4,d10.2,3d12.4,d9.1)
      end
      subroutine diffun (n, t, y, ydot)
      implicit real*8 (a-h,o-z)
      dimension y(3), ydot(3)
c
      ydot(1) = -0.04d0 * y(1) + 1.0d4 * y(2) * y(3)
      ydot(3) = 3.0d7 * y(2)**2
      ydot(2) = -ydot(1) - ydot(3)
      return
      end
      subroutine jacob (n, t, y, pd, n0)
      implicit real*8 (a-h,o-z)
      dimension y(3), pd(n0,3)
c
      pd(1,1) = -0.04d0
      pd(1,2) = 1.0d4 * y(3)
      pd(1,3) = 1.0d4 * y(2)
      pd(2,1) = 0.04d0
      pd(2,3) = -pd(1,3)
      pd(3,1) = 0.0d0
      pd(3,2) = 6.0d7 * y(2)
      pd(3,3) = 0.0d0
      pd(2,2) = -pd(1,2) - pd(3,2)
      return
      end