Most biochemical processes involve macromolecules in solution. The corresponding electrostatics is of central importance for understanding their structures and functions. An accurate and efficient numerical scheme is introduced to evaluate the corresponding electrostatic potential and force by solving the governing Poisson-Boltzmann equation. This paper focuses on the following issues: (i) the point charge singularity problem, (ii) the dielectric discontinuity problem across a molecular surface, and (iii) the infinite domain problem. Green's function associated with the point charges plus a harmonic function is introduced as the zeroth order approximation to the solution to solve the point charge singularity problem. A jump condition capturing finite difference scheme is adopted to solve the discontinuity problem across molecule surfaces, where a body-fitting grid is used. The infinite domain problem is solved by mapping the outer infinite domain into a finite domain. The corresponding stiffness matrix is symmetric and positive definite, therefore, fast algorithm such as preconditioned conjugate gradient method can be applied for inner iteration. Finally, the resulting scheme is second order accurate for both the potential and its gradient.