      PROGRAM Poisson
!Solves iteratively the Equation 
! d^2 psi / d x^2  =  S (x)
!with the boundary condition psi(0) = psi(1) = 0
!which comes from the minimization of the functional
! E = int dx [ ( psi'(x) )^2 / 2 - S(x) psi(x) ]
      IMPLICIT NONE
      INTEGER           :: i
      INTEGER           :: j
      INTEGER           :: M = 10000!No. of iterations
      INTEGER,PARAMETER :: N=100!discretization
      REAL*16           :: psi(0:N+1),S(0:N+1)
      REAL*16           :: h
      REAL*16           :: aux1
      REAL*16           :: aux2
      REAL*16           :: aux3
      OPEN(100,file='difference.dat',status='unknown')
      psi(0)   = 0.d0!Boundary conditions
      psi(N+1) = 0.d0!Boundary conditions
      h        = 1.q0/real(N+1,16)!dx
      do j=0,N+1
         S(j) = h*h * real(j,16)*h * 6.q0 !here, S(x) = 6*x. Notice we are actually defining S*h^2 for later convenience
      enddo
      do j=1,N
         psi(j) = 0.q0*real(j,16)*h*(1.q0-real(j,16)*h)!initial condition
      enddo
      do i=1,M
         do j=1,N
            psi(j) = ( psi(j+1) + psi(j-1) + S(j) ) / 2.q0 !update psi
         enddo
         aux1 = 0.q0
         do j=1,N+1
            aux1 = aux1 +  (psi(j) - psi(j-1))**2.q0 !compute kinetic energy
         enddo
         aux1 = aux1 / (2.q0 * h)
         aux2 = 0.q0
         do j=1,N+1
            aux2 = aux2 + S(j)*psi(j) !compute potential energy
         enddo
         aux2 = - aux2 / h
         aux1 = aux1 + aux2!This is the value of the energy functional
         aux3 = 0.q0
         do j=1,N
            aux2 = h*real(j,16)!this is = x
            aux3 = aux3 + abs( psi(j) - aux2*(1.q0 - aux2*aux2) )
         enddo
         aux3 = aux3*h!this is the integrated field difference from the theoretical value
         write(100,*)i,abs(aux1+0.4q0),aux3 !absolute difference from the theoretical value
      enddo
      CLOSE(100)
      OPEN(100,file='psi.dat',status='unknown')
      do i=0,N+1
         write(100,*)real(i,16)*h,psi(i)
      enddo
      CLOSE(100)
      STOP
      END PROGRAM Poisson
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!ENDO OF PROGRAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
