I am trying to replicate results of Grigoriu's Solution Of Boundary Value Problem By Monte Carlo Simulation. Essentially, we are looking for the local solution of a PDE and for this, we take samples of walkers starting at the point for which we are seeking the solution. The Boundaries are of Dirichlet and Neumann type. Walkers exit the domain only at Dirichlet boundaries and are reflected at Neumann boundaries.
Now here is where my question begins. a Brownian motion reflected at zero is
$$ X(t)=|B(t)|=\hat{B}(t)+L(t) $$
where
$$ L(t)=\lim _{\varepsilon \downarrow 0} \frac{1}{2 \varepsilon} \int_{0}^{t} 1_{(-\varepsilon, \varepsilon)}(B(s)) d s $$
is local time process. Now here I have hard time to understand the local time in order to implement it. Our walker is defined as
$$
\begin{array}{l}{X_{1}=\left|B_{1}\right|=\hat{B}_{1}+L_{1}} \\ {X_{2}=\left|B_{2}\right|=\sqrt{3}\left(\hat{B}_{2}+L_{2}\right)}\end{array}
$$
Looking online, this process gives the amount of time that the walker spends close to a Neumann boundary. I tried to add dt to account for time whenever walker crosses a Neumann boundary to reflect it back in, but it was not reflecting it back. Also, some texts were counting the number of times that walker crosses the boundary as local time, which kicks the walker outside of the domain all together. I would appreciate any comments or clarifications that helps to implement it. Thanks!