c===========================================================
c     Program demonstrating use of 'lsoda' to evaluate 
c     a definite integral. 
c
c     Also demonstrates use of optional inputs, in this 
c     case the maximum number of internally defined steps
c     allowed during one call to the solver.
c
c     usage: integral <xs> <xf> [<tol>]
c===========================================================
      program        integral

      implicit       none

      integer        iargc
      real*8         r8arg

      real*8         r8_never
      parameter    ( r8_never = -1.0d-60 )

c-----------------------------------------------------------
c     Command line arguments: integration limits and LSODA
c     (absolute) error tolerance---use a stringent default 
c     tolerance.
c-----------------------------------------------------------
      real*8         xs,         xf,       tol

      real*8         default_tol
      parameter    ( default_tol = 1.0d-12 )

c-----------------------------------------------------------
c     LSODA Variables.
c-----------------------------------------------------------
      external       fcn,        jac

      integer        neq
      parameter    ( neq = 1 )

      real*8         y(neq)
      integer        itol
      real*8         rtol,       atol
      integer        itask,      istate,     iopt
      integer        lrw

      parameter    ( lrw = 22 + neq * 16 )
      real*8         rwork(lrw)

      integer        liw
      parameter    ( liw = 20 + neq )
      integer        iwork(liw)
      integer        jt

c-----------------------------------------------------------
c     Note: Default value for 'mxstep' ('iwork(6)') is 500.
c-----------------------------------------------------------
      integer        mxstep
      parameter    ( mxstep = 50 000 )

      integer        i

c-----------------------------------------------------------
c     Parse command line arguments (initial values) ...
c-----------------------------------------------------------
      if( iargc() .lt. 2 ) go to 900

      xs   = r8arg(1,r8_never)
      if( xs .eq. r8_never ) go to 900
      xf   = r8arg(2,r8_never)
      if( xf .eq. r8_never ) go to 900
      tol  = r8arg(3,default_tol)

c-----------------------------------------------------------
c     Use pure absolute control.
c-----------------------------------------------------------
      itol  = 1
      rtol  = 0.0d0
      atol  = tol

      itask = 1 
c-----------------------------------------------------------
c     Set the optional inputs as well as the flag which 
c     tells LSODA optional inputs are being used.  A value
c     of 0 or 0.0d0 for any of the optional inputs tells
c     LSODA to use the internal default.
c-----------------------------------------------------------
      do i = 5 , 10
        iwork(i) = 0
        rwork(i) = 0.0d0
      end do
      iwork(6) = mxstep
      iopt  = 1

c-----------------------------------------------------------
c     Have LSODA compute the Jacobian numerically if 
c     necessary (it won't be in this case!)
c-----------------------------------------------------------
      jt    = 2

c-----------------------------------------------------------
c     Initialize the integral. 
c-----------------------------------------------------------
      y(1) = 0.0d0

c-----------------------------------------------------------
c     Integrate from x = xs to x = xf.  Note that LSODA
c     overwrites 'xs' with x-value in use at end of 
c     integration (normally 'xf'). 
c-----------------------------------------------------------
      istate = 1

      call lsoda(fcn,neq,y,xs,xf,
     &           itol,rtol,atol,itask,    
     &           istate,iopt,rwork,lrw,iwork,liw,jac,jt)

c-----------------------------------------------------------
c     Check return code, write result to standard output if 
c     integration was successful, or message to standard 
c     error otherwise.
c-----------------------------------------------------------
      if( istate .ge. 0 ) then
         write(*,*) y(1)
      else 
         write(0,*) 'integral:  Error return ', istate,
     &              ' from LSODA'
      end if

c-----------------------------------------------------------
c     Normal exit.
c-----------------------------------------------------------
      stop

c-----------------------------------------------------------
c     Usage exit.
c-----------------------------------------------------------
 900  continue
         write(0,*) 'usage: integral <xs> <xf> [<tol>]'
      stop
      end