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[ x^1_i, x^2_i, x^3_i, \dots x^N_i \right]$ 

  - 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}$