subroutine poisson(n,h,q,den,rv,a,b,ecoul)
c  use Numerov algorithm to solve poisson equation
c  den(n) is electron density * (4*pi*r**2)
c  rv(n) is returned as electrostatic potential * r
c  q is total electron charge
c  a(n) and b(n) are work arrays
c  ecoul is the coulomb interation energy
       implicit real*8 (a-h,o-z)
       dimension den(n),rv(n),a(n),b(n)
       sd=2
       sl=-1
       do i=1,n
       a(i)=h*den(i)/(6*i)
       enddo
       rv(1)=10*a(1)+a(2)
       do i=2,n-1
       rv(i)=10*a(i)+a(i+1)+a(i-1)
       enddo
       rv(n)=10*a(n)+a(n-1)+2*q
c use Thomas's algorithm for inverting matrix
c   Dale U. von Rosenberg, "Methods for the Numerical Solution of 
c     Partial Differential Equations, Am. Elsevier Pub., 1969, pg. 113
      a(1)=sd
      ss2=sl*sl
      do i=2,n
      a(i)=sd-ss2/a(i-1)
      enddo
      b(1)=rv(1)/sd
      do i=2,n
      b(i)=(rv(i)-sl*b(i-1))/a(i)
      enddo
      rv(n)=b(n)
      do i=n-1,1,-1
      rv(i)=b(i)-sl*rv(i+1)/a(i)
      enddo
c
c
c  calculate ecoul
c
       do i=1,n
       a(i)=rv(i)/i
       enddo
       th=1
c       ecoul=0.5d0*ddot(n,den,1,a,1)
       ecoul=0.5d0*overlap(n,th,den,a)
       write(6,*) 'ecoul = ',ecoul
c
      return
      end