  We review two common numerical schemes for Coulomb potential evaluation that
differ only in their radial part of the solutions in the spherical harmonic
expansion (SHE). One is based on finite-difference method (FDM) while the other
is based on the Green's function (GF) solution to the radial part of the
Poisson equation. We analyze the methods and observe that the FDM-based
approach appears to be more efficient in terms of the convergence with the
number of radial points, particularly for monopole (l=0). However, as a known
issue, it suffers from error accumulation as the system size increases. We
identify the source of error that comes mainly from l=1 (and sometimes l=2)
contribution of SHE induced by the charge partitioning. We then propose a
hybrid scheme by combining the two methods, where the radial solution for l=0
is obtained using the FDM method and treating the remaining terms using GF
approach. The proposed hybrid method is subsequently applied to a variety of
systems to examine its performance. The results show improved accuracy than
earlier numerical schemes in all cases. We also show that, even with a generic
set of radial grid parameters, accurate energy differences can be obtained
using a numerical Coulomb solver in standard density functional studies. ~
