program sqwell

      implicit none

      INTERFACE
         FUNCTION f(x)
            REAL f
            REAL, INTENT(in) :: x
         END FUNCTION f
      END INTERFACE

      integer, parameter :: n=500
      integer :: i,j,INFO, LWORK = 3*n-1
      real , dimension(n,n) :: A
      real , dimension(n) :: B
      real :: x, h, h2, p 
      integer , dimension(3*n-1) :: WORK

      h = 1./real(n+1)
      h2 = h*h

      A = 0.
      x = h
      A(1,1) = 2. + f(x)*h*h
      A(1,2) = -1.
      A(n,n-1) = -1.
      x = 1.-h
      A(n,n) = 2 + f(x)*h*h

      do i=2,n-1
         x = real(i)*h
         A(i,i-1) = -1.
         A(i,i) = 2. + f(x)*h*h
         A(i,i+1) = -1.
      enddo


      call SSYEV('V','U',n,A,n,B,WORK,LWORK,INFO)


      Print *, "Solved the system.  INFO = ", INFO
      Print *, "INFO = 0 is success, INFO > 0 means failure to converge"

      write(100,*)"# mode  eigenvalue"
      do i=1,n
        write(100,*)i, B(i)/(h*h)
      enddo

      do j=1,10
         write(100+j,*)"# position psi psi^2 for mode",j
         write(100+j,*) 0, 0, 0
         do i = 1,n
            x = real(i)*h
            p = A(i,j)
            write(100+j,*) x, p, p*p
         enddo
         write(100+j,*) 1, 0, 0 
      enddo

      end program sqwell

      FUNCTION f(x)
      implicit none
      real , intent(in) :: x
      real :: f

      f=0.

      END FUNCTION f