In this paper, we present a numerical scheme based on a collocation method to solve stochastic non-linear Poisson–Boltzmann equations (PBE). This equation is a generalized version of the non-linear Poisson–Boltzmann equations arising from a form of biomolecular modeling to the stochastic case. Applying the collocation method based on radial basis functions (RBFs) allows us to deal with the difficulties arising from the complexity of the domain. To indicate the accuracy of the RBF method, we present numerical results for two-dimensional models, we also study the stability of this method numerically. We examine our results with the RBF-reference value and the Chebyshev Spectral Collocation (CSC) method. Furthermore, we discuss finding the appropriate shape parameter to obtain an accurate numerical solution besides greatest stability. We have exerted the Newton–Raphson approach for solving the system of non-linear equations resulting from discretization by the RBF technique.