Background
I want to approximate an SPDE of adensity process $V_t$. The Brownian motion part is given (up to a certain precision level $\epsilon>0$) by a brownian trajectory $t\mapsto w_t$. The starting point being $V_0(x)=\delta_{x_0}(x)$ (i.e. $V_0(x)=1$ if $x=x_0$) and $V^{(n)}=u_{n\epsilon}$, (considering that we have a Brownian path) I'm given the following deterministic PDE
$$ du_t(x)=-\mu(t,x,L^{(n-1)})\partial_x u_t(x)dt + \frac{1}{2}\sigma(t,x)\rho(t,L^{(n-1)})\partial_{xx}u_t(x)dt\\ -\sigma(t,x)\sqrt{1-\rho(t,L^{(n-1)})^2}\partial_x u_t(x)dw_t, $$ and $L^{(n)}=1-\int_0^\infty V^{(n)}(x)dx, for $x>0$.
Problem
I'm in need of an approach to simulate this in R. I already looked up the package Sim.DiffProc and also others like simecol or the general packages for solving PDEs and ODEs didn't really help me (I don't have that much experience in simulating this).
Question I kindly ask you to either provide me with information on how to approach solving the (determinisitc) PDE in R or with some advice on which package I could use (considering that the parameters $\sigma, \mu, \rho$ are dependent on $L^{(n-1)}$.
Thank you in advance.
P.S. Although I know that this might belong to other communities, I'm specifically interested in the numerical solution/mathematical approach of this thing so I supposed it was right on MSE.