c-----------------------------------------------------------------------
c demonstration program for the lsodar package.
c this is the january 27, 1982 version.
c this version is in double precision.
c
c for computer systems requiring a program card, the following (with
c the c in column 1 removed) may be used..
c     program lardem(larout,tape6=larout)
c
c the lsodar package is used to solve two simple problems,
c one nonstiff and one intermittently stiff.
c if the errors are too large, or other difficulty occurs,
c a warning message is printed.  all output is on unit lout = 6.
c-----------------------------------------------------------------------
      external f1, gr1, f2, jac2, gr2
      integer iopt, iout, istate, itask, itol, iwork, jroot, jt,
     1   kroot, leniw, lenrw, liw, lrw, lout, neq, nerr, ng,
     2   nfe, nfea, nge, nje, nst
      double precision atol, er, ero, errt, rtol, rwork,
     1   t, tout, tzero, y, yt
      dimension y(2), atol(2), rwork(57), iwork(22), jroot(2)
      data lout/6/
c
      nerr = 0
c-----------------------------------------------------------------------
c first problem.
c the initial value problem is..
c   dy/dt = ((2*log(y) + 8)/t - 5)*y,  y(1) = 1,  1 .le. t .le. 6
c the solution is  y(t) = exp(-t**2 + 5*t - 4)
c the two root functions are..
c   g1 = ((2*log(y)+8)/t - 5)*y (= dy/dt)  (with root at t = 2.5),
c   g2 = log(y) - 2.2491  (with roots at t = 2.47 and 2.53)
c-----------------------------------------------------------------------
c set all input parameters and print heading.
      neq = 1
      y(1) = 1.0d0
      t = 1.0d0
      tout = 2.0d0
      itol = 1
      rtol = 0.0d0
      atol(1) = 1.0d-6
      itask = 1
      istate = 1
      iopt = 0
      lrw = 44
      liw = 21
      jt = 2
      ng = 2
      write (lout,110) itol,rtol,atol(1),jt
 110  format(1h1/1x,41h demonstration program for lsodar package/////
     1  1x,14h first problem///
     2  1x,54h problem is  dy/dt = ((2*log(y)+8)/t - 5)*y,  y(1) = 1//
     3  1x,41h solution is  y(t) = exp(-t**2 + 5*t - 4)//
     4  1x,21h root functions are../
     5  10x,30h g1 = dy/dt  (root at t = 2.5)/
     6  10x,55h g2 = log(y) - 6.2491  (roots at t = 2.47 and t = 2.53)//
     7  1x,7h itol =,i3,9h   rtol =,e10.1,9h   atol =,e10.1/
     8  1x,5h jt =,i3/////)
c
c call lsodar in loop over tout values 2,3,4,5,6.
      ero = 0.0d0
      do 180 iout = 1,5
 120    continue
        call lsodar(f1,neq,y,t,tout,itol,rtol,atol,itask,istate,
     1     iopt,rwork,lrw,iwork,liw,jdum,jt,gr1,ng,jroot)
c
c print y and error in y, and print warning if error too large.
        yt = dexp(-t*t + 5.0d0*t - 4.0d0)
        er = y(1) - yt
        write (lout,130) t,y(1),er
 130    format(1x,7h at t =,e15.7,5x,3hy =,e15.7,5x,7herror =,e12.4)
        if (istate .lt. 0) go to 185
        er = dabs(er)/atol(1)
        ero = dmax1(ero,er)
        if (er .lt. 1000.0d0) go to 140
        write (lout,135)
 135    format(//1x,41h warning.. error exceeds 1000 * tolerance//)
        nerr = nerr + 1
 140    continue
        if (istate .ne. 3) go to 175
c
c if a root was found, write results and check root location.
c then reset istate to 2 and return to lsodar call.
        write (lout,150) t,jroot(1),jroot(2)
 150    format(/1x,18h root found at t =,e15.7,5x,7hjroot =,2i5)
        if (jroot(1) .eq. 1) errt = t - 2.5d0
        if (jroot(2) .eq. 1 .and. t .le. 2.5d0) errt = t - 2.47d0
        if (jroot(2) .eq. 1 .and. t .gt. 2.5d0) errt = t - 2.53d0
        write (lout,160) errt
 160    format(1x,31h error in t location of root is,e12.4//)
        if (dabs(errt) .lt. 1.0d-3) go to 170
        write (lout,165)
 165    format(//1x,36h warning.. root error exceeds 1.0e-3//)
        nerr = nerr + 1
 170    istate = 2
        go to 120
c
c if no root found, increment tout and loop back.
 175    tout = tout + 1.0d0
 180    continue
c
c problem complete.  print final statistics.
 185  continue
      if (istate .lt. 0) nerr = nerr + 1
      nst = iwork(11)
      nfe = iwork(12)
      nje = iwork(13)
      nge = iwork(10)
      lenrw = iwork(17)
      leniw = iwork(18)
      nfea = nfe
      if (jt .eq. 2) nfea = nfe - neq*nje
      write (lout,190) lenrw,leniw,nst,nfe,nfea,nje,nge,ero
 190  format(//1x,32h final statistics for this run../
     1  1x,13h rwork size =,i4,15h   iwork size =,i4/
     2  1x,18h number of steps =,i5/
     3  1x,18h number of f-s   =,i5/
     4  1x,18h (excluding j-s) =,i5/
     5  1x,18h number of j-s   =,i5/
     6  1x,18h number of g-s   =,i5/
     7  1x,16h error overrun =,e10.2)
c
c-----------------------------------------------------------------------
c second problem (van der pol oscillator).
c the initial value problem is..
c   dy1/dt = y2,  dy2/dt = 100*(1 - y1**2)*y2 - y1,
c   y1(0) = 2,  y2(0) = 0,  0 .le. t .le. 200
c the root function is  g = y1.
c an analytic solution is not known, but the zeros of y1 are known
c to 15 figures for purposes of checking the accuracy.
c-----------------------------------------------------------------------
c set tolerance parameters and print heading.
      itol = 2
      rtol = 1.0d-6
      atol(1) = 1.0d-6
      atol(2) = 1.0d-4
      write (lout,200) itol,rtol,atol(1),atol(2)
 200  format(1h1/1x,40h second problem (van der pol oscillator)///
     1  1x,56h problem is dy1/dt = y2,  dy2/dt = 100*(1-y1**2)*y2 - y1/
     2  1x,33h            y1(0) = 2,  y2(0) = 0//
     3  1x,25h root function is  g = y1//
     4  1x,7h itol =,i3,9h   rtol =,e10.1,9h   atol =,2e10.1///)
c
c loop over jt = 1, 2.  set remaining parameters and print jt.
      do 290 jt = 1,2
      neq = 2
      y(1) = 2.0d0
      y(2) = 0.0d0
      t = 0.0d0
      tout = 20.0d0
      itask = 1
      istate = 1
      iopt = 0
      lrw = 57
      liw = 22
      ng = 1
      write (lout,210) jt
 210  format(////////1x,19h solution with jt =,i2////)
c
c call lsodar in loop over tout values 20,40,...,200.
      do 270 iout = 1,10
 220    continue
        call lsodar(f2,neq,y,t,tout,itol,rtol,atol,itask,istate,
     1     iopt,rwork,lrw,iwork,liw,jac2,jt,gr2,ng,jroot)
c
c print y1 and y2.
        write (lout,230) t,y(1),y(2)
 230    format(1x,7h at t =,e15.7,5x,4hy1 =,e15.7,5x,4hy2 =,e15.7)
        if (istate .lt. 0) go to 275
        if (istate .ne. 3) go to 265
c
c if a root was found, write results and check root location.
c then reset istate to 2 and return to lsodar call.
        write (lout,240) t
 240    format(/1x,18h root found at t =,e15.7)
        kroot = ifix(sngl(t/81.2d0 + 0.5d0))
        tzero = 81.17237787055d0 + dfloat(kroot-1)*81.41853556212d0
        errt = t - tzero
        write (lout,250) errt
 250    format(1x,31h error in t location of root is,e12.4//)
        if (errt .lt. 1.0d0) go to 260
        write (lout,255)
 255    format(//1x,33h warning.. root error exceeds 1.0//)
        nerr = nerr + 1
 260    istate = 2
        go to 220
c
c if no root found, increment tout and loop back.
 265    tout = tout + 20.0d0
 270    continue
c
c problem complete.  print final statistics.
 275  continue
      if (istate .lt. 0) nerr = nerr + 1
      nst = iwork(11)
      nfe = iwork(12)
      nje = iwork(13)
      nge = iwork(10)
      lenrw = iwork(17)
      leniw = iwork(18)
      nfea = nfe
      if (jt .eq. 2) nfea = nfe - neq*nje
      write (lout,280) lenrw,leniw,nst,nfe,nfea,nje,nge
 280  format(//1x,32h final statistics for this run../
     1  1x,13h rwork size =,i4,15h   iwork size =,i4/
     2  1x,18h number of steps =,i5/
     3  1x,18h number of f-s   =,i5/
     4  1x,18h (excluding j-s) =,i5/
     5  1x,18h number of j-s   =,i5/
     6  1x,18h number of g-s   =,i5)
 290  continue
c
c
      write (lout,300) nerr
 300  format(////1x,31h number of errors encountered =,i3)
      stop
      end
      subroutine f1 (neq, t, y, ydot)
      integer neq
      double precision t, y, ydot
      dimension y(1), ydot(1)
      ydot(1) = ((2.0d0*dlog(y(1)) + 8.0d0)/t - 5.0d0)*y(1)
      return
      end
      subroutine gr1 (neq, t, y, ng, groot)
      integer neq, ng
      double precision t, y, groot
      dimension y(1), groot(2)
      groot(1) = ((2.0d0*dlog(y(1)) + 8.0d0)/t - 5.0d0)*y(1)
      groot(2) = dlog(y(1)) - 2.2491d0
      return
      end
      subroutine f2 (neq, t, y, ydot)
      integer neq
      double precision t, y, ydot
      dimension y(2), ydot(2)
      ydot(1) = y(2)
      ydot(2) = 100.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1)
      return
      end
      subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
      integer neq, ml, mu, nrowpd
      double precision t, y, pd
      dimension y(2), pd(nrowpd,2)
      pd(1,1) = 0.0d0
      pd(1,2) = 1.0d0
      pd(2,1) = -200.0d0*y(1)*y(2) - 1.0d0
      pd(2,2) = 100.0d0*(1.0d0 - y(1)*y(1))
      return
      end
      subroutine gr2 (neq, t, y, ng, groot)
      integer neq, ng
      double precision t, y, groot
      dimension y(2), groot(1)
      groot(1) = y(1)
      return
      end