The Poisson-Boltzmann Equation (PBE) is a widely used implicit solvent model for the electrostatic analysis of solvated biomolecules. The numerical solution of the PBE is known to be challenging, due to the consideration of discontinuous coefficients, complex geometry of protein structures, singular source terms, and strong nonlinearity. In this talk, I will offer a brief overview of recent developments in the literature for resolving these numerical difficulties. For treating dielectric interface and complex geometry, both finite element methods and Cartesian grid finite difference methods have been developed for delivering a second order accuracy in space. In the framework of pseudo-time integration, we have constructed an analytical treatment to suppress the nonlinear instability. Recently, we have introduced a new regularization method for treating charge singularity in solvated biomolecules. In a regularization method, by decomposing the potential function into two or three components, the singular component can be analytically represented by the Green's function, while other components possess a higher regularity. Our new regularization combines the efficiency of two-component schemes with the accuracy of the three-component schemes. Numerical experiments of several benchmark examples and free energy calculations of protein systems are conducted to demonstrate the stability, accuracy, and efficiency of the new algorithms.