In this paper we study probabilistic and neural network approximations for solutions to Poisson equation subject to H\" older or $C^2$ data in general bounded domains of $\mathbb{R}^d$. We aim at two fundamental goals. The first, and the most important, we show that the solution to Poisson equation can be numerically approximated in the sup-norm by Monte Carlo methods based on a slight change of the walk on spheres algorithm. This provides estimates which are efficient with respect to the prescribed approximation error and without the curse of dimensionality. In addition, the overall number of samples does not not depend on the point at which the approximation is performed. As a second goal, we show that the obtained Monte Carlo solver renders ReLU deep neural network (DNN) solutions to Poisson problem, whose sizes depend at most polynomially in the dimension $d$ and in the desired error. In fact we show that the random DNN provides with high probability a small approximation error and low polynomial complexity in the dimension.