As far as literature searches go in the direction of cataloguing finite difference methods explicitly, there appears to be no descriptive representation of reaction–diffusion equations in their full three spatial dimensions. It is therefore our task (albeit tedious) to assemble the matrices and vectors for the model.

We extend the two-dimensional operator formulation to the three-dimensional case by introducing an additional spatial index \( n \in \{0,1,\dots,L\} \) with step size \( \Delta z \). The flux operators \( (\epsilon^{+}_x + \epsilon^{+}_y) \) are augmented by the \( z \)-direction term \( \epsilon^{+}_z \), so that the discrete diffusion operator acting on the substrate \( S \) is given by the 3D seven-point stencil.

3D Stencil Formulation

\[ \begin{aligned} (\epsilon^{+}_x + \epsilon^{+}_y + \epsilon^{+}_z) S^k_{l,m,n} &= a^{-}_{l,m,n;x}\, S^{k+1}_{l-1,m,n} + a^{+}_{l,m,n;x}\, S^{k+1}_{l+1,m,n} \\ &\quad + a^{-}_{l,m,n;y}\, S^{k+1}_{l,m-1,n} + a^{+}_{l,m,n;y}\, S^{k+1}_{l,m+1,n} \\ &\quad + a^{-}_{l,m,n;z}\, S^{k+1}_{l,m,n-1} + a^{+}_{l,m,n;z}\, S^{k+1}_{l,m,n+1} \\ &\quad + b_{l,m,n}\, S^{k+1}_{l,m,n}. \end{aligned} \tag{1} \]

Coefficient Definitions

\[ \begin{aligned} a^{-}_{l,m,n;x} &= \frac{\Delta t}{(\Delta x)^2}\, D\!\left(\frac{M^k_{l-1,m,n} + M^k_{l,m,n}}{2}\right),\\ a^{+}_{l,m,n;x} &= \frac{\Delta t}{(\Delta x)^2}\, D\!\left(\frac{M^k_{l,m,n} + M^k_{l+1,m,n}}{2}\right),\\[6pt] a^{-}_{l,m,n;y} &= \frac{\Delta t}{(\Delta y)^2}\, D\!\left(\frac{M^k_{l,m-1,n} + M^k_{l,m,n}}{2}\right),\\ a^{+}_{l,m,n;y} &= \frac{\Delta t}{(\Delta y)^2}\, D\!\left(\frac{M^k_{l,m,n} + M^k_{l,m+1,n}}{2}\right),\\[6pt] a^{-}_{l,m,n;z} &= \frac{\Delta t}{(\Delta z)^2}\, D\!\left(\frac{M^k_{l,m,n-1} + M^k_{l,m,n}}{2}\right),\\ a^{+}_{l,m,n;z} &= \frac{\Delta t}{(\Delta z)^2}\, D\!\left(\frac{M^k_{l,m,n} + M^k_{l,m,n+1}}{2}\right). \end{aligned} \] \[ b_{l,m,n} = 1 – \left( a^{-}_{l,m,n;x} + a^{+}_{l,m,n;x} + a^{-}_{l,m,n;y} + a^{+}_{l,m,n;y} + a^{-}_{l,m,n;z} + a^{+}_{l,m,n;z} \right). \tag{2} \]

3D Block-Matrix Assembly

\[ \mathcal{A}^k_Y = \begin{bmatrix} A^k_Y & H^{+}_1 & 0 & \cdots & 0 \\ H^{-}_2 & A^k_Y & H^{+}_2 & \cdots & 0 \\ 0 & H^{-}_3 & A^k_Y & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & H^{+}_{L-1} \\ 0 & \cdots & 0 & H^{-}_L & A^k_Y \end{bmatrix}. \tag{3} \]

Global Coupled System

\[ \mathcal{M}^k = \begin{bmatrix} \mathcal{A}^k_S & \mathcal{Z}^k_1 & 0 & \mathcal{H}^k \\ \mathcal{G}^k_1 & \mathcal{A}^k_X & 0 & 0 \\ 0 & \mathcal{Z}^k_2 & \mathcal{A}^k_I & 0 \\ \mathcal{G}^k_2 & 0 & 0 & \mathcal{A}^k_E \end{bmatrix}. \tag{4} \]

Kronecker Representation

\[ \mathcal{A}^k_Y = I_z \otimes A^k_Y + T_z \otimes I_{xy}, \tag{5} \] where \( I_z \) is the identity in the \(z\)-direction, \( I_{xy} \) is the identity in the \((x,y)\)-plane, and \( T_z \) contains the coefficients \( a^{\pm}_{l,m,n;z} \). This form preserves sparsity and allows direct use of sparse BLAS.