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