program poisson

      implicit none

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

      integer, parameter :: n=100
      integer :: i,INFO
      real , dimension(n,n) :: A
      real , dimension(n) :: B
      real :: x, h, h2, V0, VN1
      integer , dimension(n) :: IPIV

      A = 0. 
      A(1,1) = -2.
      A(1,2) = 1.
      A(n,n-1) = 1.
      A(n,n) = -2

      do i=2,n-1
         A(i,i-1) = 1.
         A(i,i) = -2.
         A(i,i+1) = 1.
      enddo

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

      V0 = 0.
      VN1 = 5.

      x = 0 + 1.*h
      B(1) = h2*f(x) - V0
      x = 1. - 1.*h
      B(n) = h2*f(x) - VN1

      do i = 2,n-1
         x = real(i)*h
         B(i) = h2*f(x)
      enddo


      call SGESV(n,1,A,n,IPIV,B,n,INFO)


      Print *, "Solved the system.  INFO = ", INFO
      Print *, "INFO = 0 is success, INFO > 0 means A is singular", INFO

      write(10,*)"# position potential "
      write(10,*) 0,  V0
      write(10,*) h,  B(1)
      do i = 2,n-1
         x = real(i)*h
         write(10,*) x, B(i)
      enddo
      write(10,*) 1-h, B(n) 
      write(10,*) 1, VN1

      call flush(10)
      end program poisson

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

      f=0.
!      f = exp(-1000.*(x-0.5)**4)

      END FUNCTION f