c-----------------------------------------------------------------------
c example problem.
c
c the following is a simple example problem, with the coding
c needed for its solution by lsode.  the problem is from chemical
c kinetics, and consists of the following three rate equations..
c     dy1/dt = -.04*y1 + 1.e4*y2*y3
c     dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
c     dy3/dt = 3.e7*y2**2
c on the interval from t = 0.0 to t = 4.e10, with initial conditions
c y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
c
c the following coding solves this problem with lsode, using mf = 21
c and printing results at t = .4, 4., ..., 4.e10.  it uses
c itol = 2 and atol much smaller for y2 than y1 or y3 because
c y2 has much smaller values.
c at the end of the run, statistical quantities of interest are
c printed (see optional outputs in the full description below).
c
      external fex, jex
      double precision atol, rtol, rwork, t, tout, y
      dimension y(3), atol(3), rwork(58), iwork(23)
      neq = 3
      y(1) = 1.d0
      y(2) = 0.d0
      y(3) = 0.d0
      t = 0.d0
      tout = .4d0
      itol = 2
      rtol = 1.d-4
      atol(1) = 1.d-6
      atol(2) = 1.d-10
      atol(3) = 1.d-6
      itask = 1
      istate = 1
      iopt = 0
      lrw = 58
      liw = 23
      mf = 21
      do 40 iout = 1,12
        call lsode(fex,neq,y,t,tout,itol,rtol,atol,itask,istate,
     1     iopt,rwork,lrw,iwork,liw,jex,mf)
        write(6,20)t,y(1),y(2),y(3)
  20    format(7h at t =,1pe12.4,6h   y =,1p3e14.6)
        if (istate .lt. 0) go to 80
  40    tout = tout*10.d0
      write(6,60)iwork(11),iwork(12),iwork(13)
  60  format(/12h no. steps =,i4,11h  no. f-s =,i4,11h  no. j-s =,i4)
      stop
  80  write(6,90)istate
  90  format(///22h error halt.. istate =,i3)
      stop
      end
 
      subroutine fex (neq, t, y, ydot)
      double precision t, y, ydot
      dimension y(3), ydot(3)
      ydot(1) = -.04d0*y(1) + 1.d4*y(2)*y(3)
      ydot(3) = 3.d7*y(2)*y(2)
      ydot(2) = -ydot(1) - ydot(3)
      return
      end
 
      subroutine jex (neq, t, y, ml, mu, pd, nrpd)
      double precision pd, t, y
      dimension y(3), pd(nrpd,3)
      pd(1,1) = -.04d0
      pd(1,2) = 1.d4*y(3)
      pd(1,3) = 1.d4*y(2)
      pd(2,1) = .04d0
      pd(2,3) = -pd(1,3)
      pd(3,2) = 6.d7*y(2)
      pd(2,2) = -pd(1,2) - pd(3,2)
      return
      end