program numerov1
C
C  illustration of Numerov algorithm and hazards thereof
C
C  integrals Poisson's equation
C
C  CWJ SDSU April 2005
C
C  solve 
C    d^2/dr^2 (r*phi) = 4*pi*r*rho
C

      implicit none

      integer size
      parameter(size=200)

      real rphi(0:size)  ! = r x phi
      real rrho(0:size)  ! = r x rho   source term

      real r
      real L            ! = size of "box"

      real dr
      integer n
      real pi
      real slop

      pi = 2.0*asin(1.0)


      print*,' enter size of box '
      read*,L

      dr = L/float(size)

      call filldensity(size,dr,rrho)
C.............. integrate outward`.........

      rphi(0)= 0.
      rphi(1) = 0.5*dr

      do n = 2,size
        rphi(n) = 2*rphi(n-1)-rphi(n-2)
     & - dr*dr/12.*4.*pi*(rrho(n)+10*rrho(n-1)+rrho(n-2))
        
        r = n*dr
        write(6,101),r,rphi(n),1-0.5*(r+2)*exp(-r),
     &  rphi(n)-(1-0.5*(r+2)*exp(-r))
  101  format(f6.3,2x,3f10.5)
        if(n.ne.0)then
           write(27,*)r,rphi(n)/r
           write(17,*)r,(1.-0.5*(r+2)*exp(-r))/r
        endif

      enddo

      end


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine filldensity(size,dr,rrho)

      implicit none

      integer size
      real rrho(0:size)
      real dr
      
      real r
      integer i
      real pi

      pi = 2.0*asin(1.)

      do i = 0,size
         r=dr*i
         rrho(i)= r*exp(-r)/8./pi
      enddo

      return
      end