The Poisson Boltzmann equation (PBE) is a well-established implicit solvent continuum model for the electrostatic analysis of solvated biomolecules. The numerical solution of the nonlinear PBE is still a challenge due to its strong singularity by the source terms, dielectrically distinct regions, and exponential nonlinear terms. In this dissertation, a new alternating direction implicit method (ADI) is proposed for solving the nonlinear PBE using a two-component regularization. This scheme inherits all the advantages of the two-component regularization and the pseudo-time solution of the PBE while possesses a novel approach to combine them. A modified Ghost Fluid Method (GFM) has been introduced to incorporate the nonzero jump condition into the ADI framework to construct a new GFM-ADI method. It produced better results in terms of spatial accuracy and stability compared to the existing ADI methods for PBE and it is simpler to implement by circumventing the work necessary to apply the rigorous 3D interface treatments with the regularization. Moreover, the stability of the GFM-ADI method has been significantly improved in comparing with the non-regularized ADI method, so that stable and efficient protein simulations can be carried out with a pretty large time step size. Two locally one-dimensional (LOD) methods have also been developed for the time-dependent regularized PBE, which are unconditionally stable. Finally, for numerical validation, we have evaluated the solvation free energy for a collection of 24 proteins with various sizes and the salt effect on the protein-protein binding energy of protein complexes.