  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. A more general convection-diffusion equation is also numerically
studied.
