I'm learning (or rather should be learning) how to solve the Neumann problem using the method of finite differences. The specific PDE I'm interested in is: $$\left\{ \begin{array}{r l} - \Delta u = f & \text{in} \, \, U \\ \frac{\partial u}{\partial n} = 0 & \text{in} \, \, \partial U. \end{array} \right.$$ where $U$ is a bounded $C^{1}$ domain in $\mathbb{R}^{d}$ and $f \in C^{\infty}(\bar{U})$.
If we discretize the domain by looking at $U \cap (n^{-1} \mathbb{Z}^{d})$ and rescale to $nU \cap \mathbb{Z}^{d}$, then we know how to solve the corresponding equation for Dirichlet boundary conditions using hitting times of the simple random walk.
Is a similar approach available for the Neumann boundary condition?