Diffusion processes can be easily interpreted with a probabilistic approach.
1D diffusion process - probabilistic approach. Let's divide the space and time in a discrete set of points and time intervals, with points of coordinates $x_n = n \Delta x$, and time instants $t_i = i \Delta t$.
Now, let's consider the probability of state transition from state $x_n$ at time instant $t_i$ to the states $x_k$ at time instant $t_{i+1}$, and let's define the probability transition as
$T(x_k,t_{i+1}; x_n, t_i) = \left\{ \begin{array} \\ d \qquad \qquad , x_k = x_{n-1} \\ 1-2d \qquad , x_k = x_{n} \\ d \qquad \qquad , x_k = x_{n+1} \\ 0 \qquad \qquad , \text{otherwise} \end{array} \right. $,
i.e. starting from $x_n$, the probability of being in state $x_{n}$ at the next time-step is $1-2d$, the probability of being in neighboring states $x_{n\pm1}$ is $d \ge 0$, and it's zero for all the other states.
Thus the overall probability of being in state $x_n$ at time $t_{i+1}$ is equal to
$p(x_n, t_{i+1}) = (1-2d) p(x_n, t_i) + d p(x_{n-1}, t_i) + d p(x_{n+1}, t_i)$,
that can be rearranged as
$p(x_n, t_{i+1}) - p(x_n, t_i) = d p(x_{n-1}, t_i) -2d p(x_n, t_i) + d p(x_{n+1}, t_i)$.
The left-hand side could be interpreted as a first-order discrete approximation of the time derivative (using explicit Euler method),
$p(x_n, t_{i+1}) - p(x_n, t_i) = \Delta t \dfrac{\partial p}{\partial t}(x_n, t_i) + o(\Delta t) $
and the right-hand side could be interpreted as a discrete approximation of the second-order space derivative
$p(x_{n-1}, t_{i}) - 2p(x_{n}, t_i) + p(x_{n+1}, t_i) = \Delta x^2 \dfrac{\partial^2 p}{\partial x^2}(x_n, t_i) + o(\Delta x^2) $,
and thus, we can rearrange the probability equation as
$ \Delta t \dfrac{\partial p}{\partial t}(x_n, t_i) + o(\Delta t) = d \Delta x^2 \dfrac{\partial^2 p}{\partial x^2}(x_n, t_i) + o(\Delta x^2)$,
and letting $\Delta x \rightarrow 0$, $\Delta t \rightarrow 0$, so that $d \frac{\Delta x^2}{\Delta t} = D$ finite,
$\dfrac{\partial p}{\partial t}(x, t) = D \dfrac{\partial^2 p}{\partial x^2}(x, t)$.
Montecarlo simulation.
I'd use a finite volume, or a finite difference method;
dividing the space domain in cells of size $\Delta x$, and using time-step $\Delta t$, so that the diffusivity $D$ and the probability $d$ are related by $D = \frac{d \Delta x^2}{\Delta t}$;
You can initialize a vector of:
the dimension of the number of the particles you're using to represent the concentration of your specie
collecting the id of the cell where you find the particles.
$\mathbf{u}_i = \left[ u^1_i, u^2_i, u^3_i, \dots u^N_i \right]$
so that the concentration in the cell $x_n$ is the number of particles of in each cell divided by the cell volume, and thus it's proportional to the number of elements of $\mathbf{u}$ equal to $x_i$, $\rho(x_n) \sim N(u^k = x_n)$.
apply transition probability $T(x_k, t_{i+1}, x_n, t_i)$, using a random or a pseudorandom number generator, to get $\mathbf{u}_{i+i}$.
Python implementation - Colab. Here, I'm attaching a Colab sheet with the implementation of the method. This is in Python, and with unnecessary for loops. Even though it's not the most efficient implementation, it should be a working one.
Here,
- 1000 particles sampled from uniform distribution in the range $x \in [0, 1]$
- infinite domain,
- pay attention at the Courant-Friedrichs-Leay (CFL) condition, that relates $D$, $\Delta x$, to the $\Delta t$ to get $d \le 0.5$, and thus probability $\in [0,1]$ in transition matrix $T$.