We will present a new Monte Carlo method for solving boundary value problems (BVPs) involving the Poisson–Boltzmann equation (PBE). Such BVPs arise in many situations involving the calculation of electrostatic properties of solvated large molecules. The PBE is one of the so-called implicit solvent models, and has accurately modeled electrostatics over a wide range of ionic solvent concentrations. With the new method we compare the algorithmic and computational properties of this algorithm to more commonly used, deterministic, techniques, and we present some computational results. In addition, recent computational enhancements have made our methods time competitive with deterministic techniques.