In this paper, a fully discrete solution method for a time-dependent diffusion-reaction system of equations based on second-order Strang splitting is proposed. The non-linear reaction operators couple the system while the diffusion and emission (right hand sides) terms are fully decoupled. Here we consider the case of anomalous (fractional) diffusion, defining the fractional power of the operators through their spectral decomposition. Fractional diffusion is non-local, which means that the corresponding discrete operators are represented by dense matrices. Dense matrix operations generally lead to a significant increase in computational complexity. To overcome this kind of difficulty, BURA (Best Uniform Rational Approximation) is applied. From an algorithmic point of view, BURA methods reduce work with dense matrices to a linear combination of solutions of linear systems with sparse matrices. The results presented further develop the construction and analysis in the earlier work [4], where the application of sequential (first-order) splitting is discussed. A composite algorithm that integrates Strang splitting, Crank-Nicolson and Runge-Kutta schemes in time and finite element methods and BURA in space is analyzed. Sufficient conditions for balancing the errors of different origin are obtained. An essential novelty of the presented results concerns the treatment of the non-local diffusion sub-problems in the second-order time-stepping case. Some readers may be surprised that matrix vector multiplication is computationally more difficult than solving a linear system in the case of fractional diffusion.