We discuss an efficient higher order finite difference algorithm for solving systems of 3-D reaction–diffusion equations with nonlinear reaction terms. The algorithm is fourth-order accurate in both the temporal and spatial dimensions. It requires only a regular seven-point difference stencil similar to that used in the standard second-order algorithms, such as the Crank–Nicolson algorithm. Numerical examples are presented to demonstrate the efficiency and accuracy of the new algorithm.