A problem with Poisson equation

Question

I'm computing the Hartree potentials of atoms by solving the Poisson equation and I use hydrogen atom as a test case. The Poisson equation for hydrogen atom in atomic units is given by $$\nabla^2 V_H = -4 \exp(-2 r)$$ where $r = \sqrt{x^2+y^2+z^2}$. The numerical solution of $V_H(x,y,z)$ with $z=0$ is illustrated in the following figure:

The numerical solution is computed with the Conjugate Gradient method so that the Laplacian is computed with stencil $$ \nabla^2 f(x,y,z) \approx \frac{a+b+c-6 f(x,y,z)}{h^2} $$ where $$a = f(x+h,y,z)+f(x-h,y,z),$$ $$b = f(x,y+h,z)+f(x,y-h,z),$$ and $$c = f(x,y,z+h)+f(x,y,z-h).$$

As the right-hand side of the Poisson equation is spherically symmetric the Poisson equation takes the form $$\frac{d^2 v_H}{dr^2} + \frac{2}{r} \frac{d v_H}{dr} = -4 \exp(-2 r)$$ where $v_H$ is the Hartree potential as a function of $r$. The solution of this equation is $$v_H(r) = - \frac{r + 1}{r} \exp(-2r),$$ which is illustrated in the following figure: (here $r=\sqrt{x^2+y^2+z^2}$ and $z=0$).

The Octave code computing the Laplacian is here:

function aLap = StencilLap3d( a, rStep )
  rH2 = rStep * rStep;
  vSize = size( a );
  nXSize = vSize( 1 );
  nYSize = vSize( 2 );
  nZSize = vSize( 3 );
  nXDim = ( nXSize - 1 ) / 2;
  nYDim = ( nYSize - 1 ) / 2;
  nZDim = ( nZSize - 1 ) / 2;
  assert( nXDim == round( nXDim ) );
  assert( nYDim == round( nYDim ) );
  assert( nZDim == round( nZDim ) );
  aLap = zeros( nXSize, nYSize, nZSize );
  rXP = 0;
  rXM = 0;
  rYP = 0;
  rYM = 0;
  rZP = 0;
  rZM = 0;
  for nX = (-nXDim):nXDim
    for nY = (-nYDim):nYDim
      for nZ = (-nZDim):nZDim
        if ( nX < nXDim )
          rXP = a( nXDim + 1 + nX + 1, nYDim + 1 + nY, nZDim + 1 + nZ );
        else
          rXP = 0.0;
        endif
        if ( nX > -nXDim )
          rXM = a( nXDim + 1 + nX - 1, nYDim + 1 + nY, nZDim + 1 + nZ );
        else
          rXM = 0.0;
        endif
        if ( nY < nYDim )
          rYP = a( nXDim + 1 + nX, nYDim + 1 + nY + 1, nZDim + 1 + nZ );
        else
          rYP = 0.0;
        endif
        if ( nY > -nYDim )
          rYM = a( nXDim + 1 + nX, nYDim + 1 + nY - 1, nZDim + 1 + nZ );
        else
          rYM = 0.0;
        endif
        if ( nZ < nZDim )
          rZP = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ + 1 );
        else
          rZP = 0.0;
        endif
        if ( nZ > -nZDim )
          rZM = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ - 1 );
        else
          rZM = 0.0;
        endif
        n0 = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ );
        aLap( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ ) = ...
          ( rXP + rXM + rYP + rYM + rZP + rZM - 6 * n0 ) / rH2;
      endfor
    endfor
  endfor
endfunction

Obviously the solution of the differential equation should be the same as the numerical solution, so there is something wrong here. Based on my previous computations I suspect that the error is in the differential equation. Can somebody tell what is wrong?

Answer 1

Problem solved: I just have to include the solution of the homogeneous equation into the solution of the differential equation.

The solution to the differential equation is $$V_H(r) = -\frac{K}{r} - \frac{r+1}{r} \exp(-2r).$$ By requiring that $\lim_{r \to 0} V_H(r)$ is real we get $K=-1$ and $$V_H(r) = \frac{1}{r} - \frac{r+1}{r} \exp(-2r).$$