Recently, we have introduced a pseudo-time coupled nonlinear partial differential equation (PDE) model for biomolecular solvation analysis. Based on a free energy optimization, a boundary value system is derived to couple a nonlinear Poisson-Boltzmann (NPB) equation for electrostatic potential with a generalized Laplace-Beltrami (GLB) equation defining the biomolecular surface. By introducing a pseudo-time in both processes, a more efficient coupling is achieved through the steady state solution of two nonlinear parabolic PDEs. For the GLB equation, the pseudo-transient continuation is attained via a potential driven geometric flow PDE, which defines a smooth biomolecular surface to characterize the dielectric boundary between biomolecules and the surrounding aqueous environment. The resulting smooth dielectric profile, however, introduce some instability issue in solving the time-dependent NPB equation. This motivates us to develop an operator splitting alternating direction implicit (ADI) scheme, in which the nonlinear instability is completely avoided through analytical integration. To speed up the computation of molecular surface, a new fully implicit ADI scheme is developed too for solving the geometric flow equation. Unconditional stability can be realized by both ADI schemes in solving unsteady NPB and GLB equations separately. In solving a coupled system for real biomolecules and chemical compounds, the proposed numerical schemes are found to be conditionally stable. Nevertheless, the time stability can be maintained by using very large time increments, so that the present biomolecular simulation becomes much faster.