  A simple-to-implement weak-sense numerical method to approximate reflected
stochastic differential equations (RSDEs) is proposed and analysed. It is
proved that the method has the first order of weak convergence. Together with
the Monte Carlo technique, it can be used to numerically solve linear parabolic
and elliptic PDEs with Robin boundary condition. One of the key results of this
paper is the use of the proposed method for computing ergodic limits, i.e.
expectations with respect to the invariant law of RSDEs, both inside a domain
in $\mathbb{R}^{d}$ and on its boundary. This allows to efficiently sample from
distributions with compact support. Both time-averaging and ensemble-averaging
estimators are considered and analysed. A number of extensions are considered
including a second-order weak approximation, the case of arbitrary oblique
direction of reflection, and a new adaptive weak scheme to solve a Poisson PDE
with Neumann boundary condition. The presented theoretical results are
supported by several numerical experiments.
