We introduce Monte Carlo methods to compute the solution of elliptic equations with pure Neumann boundary conditions. We first prove that the solution obtained by the stochastic representation has a zero mean value with respect to the invariant measure of the stochastic process associated to the equation. Pointwise approximations are computed by means of standard and new simulation schemes especially devised for local time approximation on the boundary of the domain. Global approximations are computed thanks to a stochastic spectral formulation taking into account the property of zero mean value of the solution. This stochastic formulation is asymptotically perfect in terms of conditioning. Numerical examples are given on the Laplace operator on a square domain with both pure Neumann and mixed Dirichlet-Neumann boundary conditions.