c=================
c   Term Project 11-21-98
c================

      program     schrodinger

      implicit    none

      
c-------
c  Mesh
c-------
      h_x  = (xmax-xmin) / 2**level
      nx  = 2**level + 1
      do i = 1, nx
         x(i) = xmin + (i-1)*h_x
      end do
c--------
c   Call Potential Function, or use default of square well
c----------
      if (vname .ne. def_vname) then
         call potential(vname, vx, x, nx, h_x)
c--------------------------------
c--- NOTE: vname : name of file containing the potential function
              vx : storage which will contain the function 
      else (((((SET UP SQUARE WELL)))))
c--------------------------------------
c   Call Initial Wavefunction or use default of a Gaussian
c-------------
      if (wfname .eq. def_wfname) then
         left_amp = 1.0d0
         del = 
         x0  =
 
         call dvgaussian (l, dl, ddl, x, nx, left_amp, x0, del
      else
         call initialize(wfname, psi, x, nx, h_x)  
c------------------------------------
c------ NOTE:  wfname : name of file containing the function. 
c                 psi : storage for the wavefunction
c---------------
      
c-----------------
c  Set up time variables
c----------------
      h_t = dt_dx * h_x
c--------------------
c   Set up initial wavefunction

c==================
c   Time Step Loop
c===================
      n    = 2
      nm1  = 1
      np1  = 1
      ofreq = 2**(level - olevel)

      call vsxynt('wave', 0.0d0, x, psi(1, np1), nx)
      call gnuout(psi(1, np1), x, nx, 0.0d0, ofreq)

      if (level .eq. olevel) then
      write(0,*) ' Olevel is equal to level'
         call vsxynt('wave', h_t, x, psi(1, np1), nx)
         call gnuout(psi(1, np1), x, nx, h_t, ofreq)
      endif
 
c-----------------
c  Temporarily
c-----------------
      nt = nx
c-----------------
c-----------------
c   Set up coefficients
c------------------
      cj1  = 1/(4.0d0 * h_x**2)

      do j = 1, nx
        cj(j)   = I/h_t - 0.50d0/h_x**2 + vx(j)/2
        crj(j)  = I/h_t + 1/(2 * h_x**2) + vx(j)/2
      end do 
c
c---NOTE:  cj1 is the coefficient of j-1 and j+1 
c   cj is coeff. of j
c   crj is coeff. of right hand side vector at j

      do j = 1, nx
        rhs(j)  = psi(j) * crj(j) + psi(j+1) * (-cj1) +
     &            psi(j-1) * (-cj1)
      end do

      call tridiagonalizer( cj1, cj, psi, rhs)
c----------
c  The point of the subroutine "tridiagonalizer is to use the coefficients
c  passed in along with the function psi(j) (for the current time step)
c  and solve the tridiagonal system that results from the new psi(j) at
c  each time step.  The point of this solver is to return psi(j,n+1) given
c  psi(j,n) and the correct coefficients.
c----------

 c------------------- c   Set up left boundary
c-------------------
      d(1)   = psi(1,0)
      du(1)  = 0.0d0
      rhs(1) = uexact(1)

c--------------------
c   Interior of tridiag system
c--------------------
    
      do j = 2, n-1
         dl(j-1) = co_jm1_n1
         d(j)    = co_j_n1(j)
         du(j)   = co_j1_n1
         rhs(j)  = psi(j)



      end !---End of Program

      subroutine  tridiagonalizer(cj1,   , psi, rhs)

      implicit    none

c--------------
c   Set up boundaries and interior
c--------------
      d(1)       = 1.0d0
      du(1)      = 0.0d0
      rhs(1)     = uexact(1)

      do j = 2, nx-1
         dl(j-1) = cj1
         d(j)    = cj(j)
         du(j+1) = cj1
         rhs(j)  = rhs(j)
      end do

      dl(nx-1)    = 0.0d0
      d(nx)       = 1.0d0
      rhs(nx)     = uexact(n)

      nrhs  = 1
      call dgtsv(nx, nrhs, dl, d, du, rhs, nx, info)

      if (info .eq. 0) then 
c Solver Successful