The neutron flux $N(\mathbf{x}, t)$ equation
\begin{align}
\dfrac{\partial{N}}{\partial{t}} &= \sigma^2 \Biggl[ \underbrace{\dfrac{\partial^2{N}}{\partial{x_1}^2} + \dfrac{\partial^2{N}}{\partial{x_2}^2} + \dfrac{\partial^2{N}}{\partial{x_3}^2}}_{\nabla^2{N} \equiv \text{“Laplacian N”}} \Biggr] \tag{1}
\\
\\
\Longrightarrow\dfrac{\partial N}{\partial t} & = \sigma^2 \nabla^2 N + A \cdot N
\end{align}
so when we setup the difference equations there are 3 time-space partial derivatives to account for whereas we had only 1 in the derivation. Each of the x, y & z dimensions of the cube of radioactive material come into play. Recalling the heat equation
$$ \underbrace{\dfrac{u(x,t + \Delta t ) \; – \; u(x,t)}{\Delta t}}_{\text{time step}} \simeq \dfrac{\partial{u}}{\partial{t}} = \dfrac{\partial^2{u}}{\partial{x^2}} \simeq \underbrace{\dfrac{u(x + \Delta x, t )\; – \; 2u(x,t) \; + \; u(x – \Delta x,t)}{(\Delta x)^2}}_{\text{time-space step in 1 dimension}}$$
However, from $(1)$ we now have a three dimensional right hand equation side. We can index $x$ in each of the three dimensions such that while
$$ h \equiv \Delta x \; , \qquad k \equiv \Delta t$$
is the same as previously, we now have this increment occuring in three dimensions and so the corresponding $x = jh$ in each dimension is $x_1 = j_1h$, $x_2 = j_2h$ and $x_3 = j_3h$ for of the three x,y,z dimensions of the cubic volume. Replacing the heat equation $u$ with $N$ and $x$ with three instances $x_1$, $x_2$ and $x_3$, the Neutron Flux equation
$$\dfrac{\partial{N}}{\partial{t}} = \dfrac{\partial^2{N}}{\partial{x_1}^2} + \dfrac{\partial^2{N}}{\partial{x_2}^2} + \dfrac{\partial^2{N}}{\partial{x_3}^2}$$
using difference equations approximates to
\begin{align}\underbrace{\dfrac{N(x,t + \Delta t ) \; – \; N(x,t)}{\Delta t}}_{\text{time step}} = \dfrac{\sigma^2}{(\Delta x)^2}
\underbrace{
\begin{cases}
N((x_1 + \Delta x_1), x_2, x_3, t )\; – \; 2N(x_1,t) \; + \; N((x_1 – \Delta x_1), x_2, x_3, t )\\
+ N(x_1 , (x_2 + \Delta x_2), x_3, t )\; – \; 2N(x_2,t) \; + \; N(x_1 , (x_2 2 \Delta x_2), x_3, t )\\
+ N(x_1 , x_2, (x_3 + \Delta x_3), t )\; – \; 2N(x_3,t) \; + \; N(x_1 , x_2, (x_3 – \Delta x_3), t )
\end{cases} \Biggr\}
}_{\text{time-space step in 3 dimensions}}
\end{align}
If we define $N_{j1}^i$ as the temperature evolution (in direction $x1$) after $i$ time steps,
\begin{align}N_{j_1}^i &\simeq N(j_1h, k), \\ \\ N_{j_2}^i &\simeq N(j_2h, k), \\ \\ N_{j_3}^i &\simeq N(j_3h, k)
\end{align}
and
$$\mathbf{N_{j}^i} = \left(\begin{array}{c}N^i_{j_1} \\ N^i_{j_2} \\ N^i_{j_3}\end{array}\right)$$
we can simplify as
\begin{align}
\dfrac{N_j^{i+1} \; – \; N^i_j}{\Delta t} = \dfrac{\sigma^2}{h^2}
\begin{cases}
N^i_{j_1 + 1, j_2, j_3} \; – \; 2\mathbf{N_{j}^i} \; + \; N^i_{j_1 – 1, j_2, j_3}\\
+ N^i_{j_1, j_2 +1 , j_3} \; – \; 2\mathbf{N_{j}^i} \; + \; N^i_{j_1, j_2-1, j_3}\\
+ N^i_{j_1, j_2, j_3 + 1} \; – \; 2\mathbf{N_{j}^i} \; + \; N^i_{j_1, j_2, j_3-1}
\end{cases} \biggr\} + A \cdot \mathbf{N_{j}^i}
\end{align}
and using more compact notation,
\begin{align}
\dfrac{N_j^{i+1} \; – \; N^i_j}{\Delta t} = \dfrac{\sigma^2}{h^2} \biggl( & N^i_{j + e_1} \; – \; 2N_{j}^i \; + \; N^i_{j – e_1}\\
&+ N^i_{j + e_2} \; – \; 2N_{j}^i \; + \; N^i_{j – e_2}\\
&+ N^i_{j + e_3} \; – \; 2N_{j}^i \; + \; N^i_{j – e_3}
\biggr) + A \cdot N_{j}^i
\end{align}
where,
$$ e_1 = \left(\begin{array}{c}1 \\ 0 \\ 0 \end{array}\right), \quad e_2 = \left(\begin{array}{c}0 \\ 2 \\ 0 \end{array}\right) \quad e_3 = \left(\begin{array}{c}0 \\ 0 \\ 1 \end{array}\right)$$