Electrostatic forces and the electrostatic properties of molecules in solution are among the most important issues in understanding the structure and function of large biomolecules. The use of implicit-solvent models, such as the Poisson-Boltzmann equation (PBE), have been used with great success as a way of computationally deriving electrostatics properties such molecules. We discuss how to solve an elliptic system of partial differential equations (PDEs) involving the Poisson and the PBEs using path-integral based probabilistic, Feynman-Kac, representations. This leads to a Monte Carlo method for the solution of this system which is specified with a stochastic process, and a score function. We use several techniques to simplify the Monte Carlo method and the stochastic process used in the simulation, such as the walk-on-spheres (WOS) algorithm, and an auxiliary sphere technique to handle internal boundary conditions. We then specify some optimizations using the error (bias) and variance to balance the CPU time. We show that our approach is as accurate as widely used deterministic codes, but has many desirable properties that these methods do not. In addition, the currently optimized codes consume comparable CPU times to the widely used deterministic codes. Thus, we have an very clear example where a Monte Carlo calculation of a low-dimensional PDE is as fast or faster than deterministic techniques at similar accuracy levels.