      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) ]
!
!In this version, iterate until a required convergence is obtained
      IMPLICIT NONE
      INTEGER           :: i
      INTEGER           :: j
      INTEGER           :: M!No. of iterations
      INTEGER,PARAMETER :: N=100!no. of points
      REAL*16           :: psi(0:N+1),S(0:N+1)
      REAL*16           :: psiold(0:N+1)
      REAL*16           :: h = 1.q0/real(N+1,16)!dx
      REAL*16           :: w,wmax,dw
      REAL*16           :: aux1
      REAL*16           :: aux2
      REAL*16           :: error
      OPEN(100,file='w-iterations-psierror.dat',status='unknown')
      psi(0)   = 0.d0!Boundary conditions
      psi(N+1) = 0.d0!Boundary conditions
      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 (when computing the psi update)
      enddo
      w    = 0.5q0
      wmax = 2.q0
      dw   = 0.05q0
      do while(w.lt.wmax)
         do j=1,N
            psi(j) = 0.q0!initial condition
         enddo
         psiold(:) = psi(:)
         error = 1.q0
         M     = 0
         do while(error.gt.1.q-6)
            do j=1,N
               aux1   = ( psi(j+1) + psi(j-1) + S(j) ) / 2.q0 
               psi(j) = (1.q0 - w)*psi(j) + w*aux1!update psi using the relaxation algorithm
            enddo
            M = M + 1
            if(mod(M,20).eq.0)then!compute error after 20 iterations for the stop criterion
               error = 0.q0
               do j=1,N
                 error = error + abs( psi(j) - psiold(j) )
               enddo
               error = error / real(N,16)
               write(6,*)M,error
               psiold(:) = psi(:)!update psiold
            endif
         enddo
         error = 0.q0
         do j=0,N+1
            aux1  = real(j,16)*h
            error = error + abs( psi(j) -  aux1*(1.q0-aux1*aux1) )
         enddo
         error = error*h
         write(100,*)w,M,error
         w = w + dw
      enddo
      CLOSE(100)
      STOP
      END PROGRAM Poisson
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!ENDO OF PROGRAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
