c===========================================================
c     tlsoda: Driver routine which uses ODEPACK routine   
c     lsoda to solve 
c
c        u''(x) = -u(x)
c===========================================================
      program         tlsoda

      implicit        none

      integer         iargc,        i4arg
      real*8          r8arg

c-----------------------------------------------------------
c     Number of ODEs (when written in canonical first order
c     form).
c-----------------------------------------------------------
      integer         neq          
      parameter     ( neq = 2 )

c-----------------------------------------------------------
c     Command-line arguments (initial values of y and 
c     dy/dx).
c-----------------------------------------------------------
      real*8          y0(neq)
      real*8          r8_never
      parameter     ( r8_never = -1.0d-60 )

c-----------------------------------------------------------
c     Vector for storing requested x-values for solution
c     output.
c-----------------------------------------------------------
      integer         mxnxout
      parameter     ( mxnxout = 10 000 )
      real*8          xout(mxnxout)
      integer         nxout,            ixout

c-----------------------------------------------------------
c     LSODA variables.   Note that 'fcn' and 'jac' are 
c     user supplied subroutines (not functions) which 
c     evaluate the RHSs of the ODEs and the Jacobian 
c     of the system.  Under normal operation, the Jacobian
c     evaluator can be a dummy routine.
c-----------------------------------------------------------
      external       fcn,        jac

      real*8         y(neq)
      real*8         xbgn,       xend

      integer        itol
c-----------------------------------------------------------
c     The following comment block is extracted from the 
c     LSODA documentation.
c-----------------------------------------------------------
c rtol   = relative tolerance parameter (scalar).
c atol   = absolute tolerance parameter (scalar or array).
c   the estimated local error in y(i) will be controlled so 
c   as to be less than
c        ewt(i) = rtol*abs(y(i)) + atol     if itol = 1, or
c        ewt(i) = rtol*abs(y(i)) + atol(i)  if itol = 2.
c   thus the local error test passes if, in each component,
c   either the absolute error is less than atol (or atol(i)),
c   or the relative error is less than rtol.
c   use rtol = 0.0 for pure absolute error control, and
c   use atol = 0.0 (or atol(i) = 0.0) for pure relative error
c   control.  CAUTION.. actual (global) errors may exceed 
c   these local tolerances, so choose them CONSERVATIVELY.
c-----------------------------------------------------------
      real*8         rtol,       atol
      integer        itask,      istate,     iopt
      integer        lrw

c-----------------------------------------------------------
c     LSODA work arrays.
c-----------------------------------------------------------
      parameter    ( lrw = 22 + neq * 16 )
      real*8         rwork(lrw)

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

c-----------------------------------------------------------
c     Error tolerances.
c-----------------------------------------------------------
      real*8         tol
      real*8         default_tol
      parameter    ( default_tol = 1.0d-6 )

c-----------------------------------------------------------
c     Locals.
c-----------------------------------------------------------
      integer         ieq

c-----------------------------------------------------------
c     Argument parsing.
c-----------------------------------------------------------
      if( iargc() .ne. 2 )  go to 900
      do ieq = 1 , neq
         y0(ieq) = r8arg(ieq,r8_never)
         if( y0(ieq) .eq. r8_never ) go to 900
      end do

c-----------------------------------------------------------
c     Get requested output x-values from standard input
c-----------------------------------------------------------
      call dvfrom('-',xout,nxout,mxnxout)

c-----------------------------------------------------------
c     Use 'default_tol' for all tolerances.
c-----------------------------------------------------------
      tol    = default_tol
c-----------------------------------------------------------
c     Set LSODA parameters ... see LSODA documentation
c     for fuller description.
c-----------------------------------------------------------

      itol   = 1         ! Indicates that 'atol' is scalar
      rtol   = tol       ! Relative tolerance
      atol   = tol       ! Absolute tolerance
      itask  = 1         ! Normal computation
      iopt   = 0         ! Indicates no optional inputs
      jt     = 2         ! Jacobian type 

c-----------------------------------------------------------
c     Initialize the solution.
c-----------------------------------------------------------
      do ieq = 1 , neq
         y(ieq) = y0(ieq)
      end do

c-----------------------------------------------------------
c     Output initial solution.
c-----------------------------------------------------------
      write(*,*) xout(1),   y(1)

c-----------------------------------------------------------
c     Loop over requested output x-values ...
c-----------------------------------------------------------
      do ixout = 2, nxout
c-----------------------------------------------------------
c        Set interval limits of integration.
c-----------------------------------------------------------
         xbgn = xout(ixout-1) 
         xend = xout(ixout)

c-----------------------------------------------------------
c        Call lsoda to integrate system on
c        x = xbgn ... xend.
c
c        'istate' should always be set to 1 prior to 
c        invocation; LSODA also uses it as a return code.
c-----------------------------------------------------------
         istate = 1

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

c-----------------------------------------------------------
c        Check return code and exit with error message if 
c        there was trouble.
c-----------------------------------------------------------
         if( istate .lt. 0 ) then
            write(0,1000) istate, ixout, nxout, xout(ixout-1),
     &                    xout(ixout)
1000        format(/' sode: Error return ',i2,
     &              ' from integrator LSODA.'/
     &              ' sode: At output time ',i5,' of ',i5/
     &              ' sode: Interval ',1pe11.3,' .. ',1pe11.3/)
            go to 500
         end if
c-----------------------------------------------------------
c        Output the solution.
c-----------------------------------------------------------
         write(*,*) xout(ixout),   y(1)
      end do

 500  continue

      stop

 900  continue
         write(0,*) 'usage: tlsoda <y_i> <dy/dx_i>'
      stop
      
      end

c===========================================================
c     Implements differential equations:
c
c     u'' = -u
c
c     y(1) := u
c     y(2) := u'
c
c     y(1)' :=  y(2)
c     y(2)' := -y(1)
c
c     Called by ODEPACK routine LSODA.
c===========================================================
      subroutine fcn(neq,t,y,yprime)
         implicit    none

         integer     neq
         real*8      t,     y(neq),    yprime(neq)

         yprime(1) =  y(2)
         yprime(2) = -y(1)

         return
      end

c===========================================================
c     Implements Jacobian (optional).  Dummy routine in 
c     this case.
c===========================================================
      subroutine jac
         implicit    none

         return
      end
