c===========================================================
c     deut: Uses LSOSA to integrate ODEs which define a 
c     simple model for a deuteron (spherically symmetric,
c     time-independent Schrodinger equation with a square
c     well potential.  
c
c     usage: deut <x0> <E> <xmax> <dxout> <tol>
c 
c         <x0>     :=  Range of the potential
c         <E>      :=  Estimate of energy eigenvalue, for 
c                      any <x0>, there is a single <E0>
c                      which results in a wave function
c                      which -> 0 as <x0> -> infinity.
c         <xmax>   :=  Maximum range of integration
c         <dxout>  :=  Output step. We use this as a 
c                      parameter rather than an output 
c                      level for ease in extending 
c                      integrations to larger <xmax> as 
c                      eigenvalue E is determined more 
c                      precisely.
c         <tol>    :=  LSODA tolerance parameter.
c
c     Output to standard output is 
c
c         <x_j>   <u(x_j)>  <du(x_j)/dx>  <sign(u(x_j))>
c
c     x_j = 0, dxout, 2 dxout, ... xmax
c    
c     See class notes and Arfken, Math. Methods for 
c     Physicists, 2nd Edition, section 9.1.2
c     for more details.
c===========================================================
      program        deut

      implicit       none

      integer        iargc
      real*8         r8arg

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

c-----------------------------------------------------------
c     Command-line arguments (Note:  x0 and E are defined in
c     fcn.inc)
c-----------------------------------------------------------
      real*8         xmax,       tol

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

      integer        neq
      parameter    ( neq = 2 )

      real*8         y(neq)
      real*8         x,          xout
      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     Common communication with routine 'fcn' in 'fcn.f'.
c-----------------------------------------------------------
      include       'fcn.inc'

c-----------------------------------------------------------
c     Locals.
c-----------------------------------------------------------
      real*8         dxout

c-----------------------------------------------------------
c     Parse command line arguments. 
c-----------------------------------------------------------
      if( iargc() .ne. 5 ) go to 900

      x0      = r8arg(1,r8_never)
      E       = r8arg(2,r8_never)
      xmax    = r8arg(3,r8_never)
      dxout   = r8arg(4,r8_never)
      tol     = r8arg(5,r8_never)
      if( x0   .eq. r8_never  .or.  E     .eq. r8_never  .or.
     &    xmax .eq. r8_never  .or.  dxout .eq. r8_never  .or.
     &    tol  .eq. r8_never ) go to 900

c-----------------------------------------------------------
c     Set LSODA parameters.  Use same value for absolute 
c     and relative tolerance.
c-----------------------------------------------------------
      itol   = 1
      rtol   = tol
      atol   = tol
      itask  = 1 
      iopt   = 0
      jt     = 2

c-----------------------------------------------------------
c     Initialize the solution, and output it.
c-----------------------------------------------------------
      x    = 0.0d0
      y(1) = 0.0d0
      y(2) = 1.0d0
      write(*,1000) x, y, int(sign(1.0d0,y(1)))
1000  format(1P,3E24.16,0p,i4)

c-----------------------------------------------------------
c     Do the integration.
c-----------------------------------------------------------
      istate = 1
      do while( x .lt. xmax )
         xout   = x + dxout
         call lsoda(fcn,neq,y,x,xout,
     &              itol,rtol,atol,itask,    
     &              istate,iopt,rwork,lrw,iwork,liw,jac,jt)

         if( istate .lt. 0 ) then
            write(0,*) 'deut:  Error return ', istate,
     &                 ' from LSODA '
            write(0,*) 'deut: Current interval ', 
     &                 x, x + dxout
            stop
         end if
c-----------------------------------------------------------
c        Output the solution.
c-----------------------------------------------------------
         write(*,1000) x, y, int(sign(1.0d0,y(1)))
      end do

      stop

900   continue
         write(0,*) 'usage: deut <x0> <E> <xmax> '//
     &              '<dxout> <tol>'
      stop

      end