c===========================================================
c     Program demonstrating use of 'lsoda' to evaluate 
c     a definite integral. 
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         xlim

      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     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     Set LSODA parameters ... use pure absolute control.
c-----------------------------------------------------------
      itol  = 1

      rtol  = 0.0d0
      atol  = tol

      itask = 1 
      iopt  = 0
      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 
