c=======================================================
c     se:  Program to solve schrodinger equation:
c     
c     -u''(x) + (V(x) - E)u = 0
c     
c     Boundary conditions?????
c
c
c     usage: se <xmin> <xmax> <u0> <du0> <tol> <olevel>
c
c
c     output to std out is
c
c     x_it   u_it   u'_it
c
c     it = 1,2,... nout
c     where:
c     x_1 = xmin
c     x_nxout = xmax
c     nxout = 2**olevel + 1
c
c=======================================================
      program se

      implicit none
c      include 'fcn.inc'

      integer iargc, i4arg
      real*8  r8arg

c-------------------------------------------------------
c     Command-line arguments:
c
c     xmin:    initial integration time
c     xmax:    final integration time
c     u0:      initial value: u(0)
c     du0:     initial value: u'(0)
c     tol:     error tolerance (this program uses LSODA's
c              pure absolute error control)
c     olevel:  output level:  
c                  dxout = (xmax-xmin)/2**olevel
c-------------------------------------------------------
      real*8      xmin,xmax,u0,du0,tol
      integer     olevel
      real*8      r8_never
      parameter ( r8_never = -1.0d-60 )


c=======================================================
c     start of LSODA declarations
c-------------------------------------------------------

c-----------------------------------------------------------
c     Note that 'fcn' and 'jac' are user supplied SUBROUTINES 
c     (not functions) which evaluate the RHSs of the ODEs and 
c     the Jacobian of the system.  Under normal operation, 
c     (as in this case), the Jacobian evaluator can be a 
c     'dummy' routine; if and when needed, LSODA will compute
c     a finite-difference approximation to the Jacobian.
c-----------------------------------------------------------
      external       fcn,        jac

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

c-----------------------------------------------------------
c     y(neq): Storage for approximate solution
c     x:      Initial time for LSODA integration sub-interval
c     xout:   Final time for LSODA integration sub-interval
c-----------------------------------------------------------
      real*8         y(neq)
      real*8         x,         xout

c-----------------------------------------------------------
c     Tolerance parameters:
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        itol

c-----------------------------------------------------------
c     Control parameters and return code (see below).
c-----------------------------------------------------------
      integer        itask,      istate,     iopt

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

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


c-----------------------------------------------------------
c     'jt' defines which type of Jacobian is supplied or
c     computed; we use jt = 2 here which, as mentioned 
c     above, instructs LSODA to compute a finite-difference
c     approximation to the Jacobian if and when needed.
c-----------------------------------------------------------
      integer        jt

c===========================================================
c     End of LSODA declarations
c===========================================================

c-----------------------------------------------------------
c     miscellaneous variables
c-----------------------------------------------------------
      real*8     dxout
      integer    it, nxout

c-----------------------------------------------------------
c     argument parsing
c-----------------------------------------------------------
      if( iargc() .ne. 6 ) go to 900
      xmin     = r8arg(1, r8_never)
      xmax     = r8arg(2, r8_never)
      u0       = r8arg(3, r8_never)
      du0      = r8arg(4, r8_never)
      tol      = r8arg(5, r8_never)
      olevel   = i4arg(6, -1)
      if (xmax .le. xmin .or. u0 .eq. r8_never .or.
     &     du0 .eq. r8_never .or. tol .eq. r8_never
     &     .or. tol .gt. 10d-2 .or. tol .lt. 10d-12)
     &     go to 800

c-----------------------------------------------------------
c     set LSODA parameters ... see LSODA documentation
c     for more complete description
c-----------------------------------------------------------
      itol    = 1      ! indicates that 'atol' is scalar
      rtol    = 0.0d0  ! use pure absolute tolerance
      atol    = tol    ! absolute tolerance
      itask   = 1      ! normal computation
      iopt    = 0      ! indicates no optional inputs
      jt      = 2      ! jacobian type

c-----------------------------------------------------------
c     compute number of outputs and output interval
c     and initialize sub-interval start time and solution
c     estimate.
c-----------------------------------------------------------
      nxout  = 2**olevel +1
      dxout  = (xmax - xmin) / (nxout -1)
      x      = xmin
      y(1)   = u0
      y(2)   = du0

c-----------------------------------------------------------
c     output initial solution
c-----------------------------------------------------------
      write(*,*) x, y(1)

c-----------------------------------------------------------
c     Loop over requested outputs ... 
c
c     Set istate to 1 to indicate initial call, istate 
c     should be set to 2 for subsequent calls, but lsoda
c     will automatically do this so long as the initial
c     call is successful.
c-----------------------------------------------------------
      istate = 1

      do it = 2, nxout
         if (x .eq. 0.0d0) then
            x = x + dxout
         end if
c     set final integration position for current interval
         xout = x + dxout


c     call lsoda to integrate system on [x ... xout]
         call lsoda(fcn, neq, y, x, xout,
     &                 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, it, nxout, x, x + dxout
1000       format(/' sode: Error return ',i2,
     &           ' from integrator LSODA.'/
     &           ' sode: At output time ',i5,' of ',i5/
     &           ' sode: Interval ',1p,e11.3,0p,
     &           ' .. ',1p,e11.3,0p/)
            go to 500
         end if
         
         write(*,*) x, y(1)
         
         
      end do

 500  continue
      stop


      stop

 800  continue
         write(0,*) 'xmax > xmin'
         write(0,*) '10e-12 <= tol <= 10e-2'
         write(0,*) 'olevel > 0'
         write(0,*)


 900  continue
         write(0,*) 'usage: se <xmin> <xmax> <u0> <du0> '//
     &              '<tol> <olevel>'
      stop
  
      end

      


      subroutine fcn(neq,x,y,yprime)
         implicit none
         
         real*8   calcv2
         external calcv2

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

         real*8  v
         
         integer     nmax
         parameter ( nmax = 100 )

         call calcv(v, x, nmax)

         
         yprime(1) = y(2)
         yprime(2) = (v - 1.0d0) * y(1)

         return
      end

      subroutine jac
         implicit none

         return
      end


