In this study, the authors present high-fidelity numerical simulations for capturing the Turing pattern in a multi-dimensional Gray–Scott reaction-diffusion system. For this purpose, an explicit mixed modal discontinuous Galerkin (DG) scheme based on multi-dimensional structured meshes is employed. This numerical scheme deals with high-order derivatives in diffusion, and highly nonlinear functions in reaction terms. Spatial discretization is accomplished using hierarchical basis functions based on scaled Legendre polynomials. A new reaction term treatment is also detailed, showing an intrinsic property of the DG scheme and preventing incorrect solutions due to excessively nonlinear reaction terms. The system is reduced to a set of time-dependent ordinary differential equations, which are solved with an explicit third-order Total Variation Diminishing (TVD) Runge–Kutta method. The Saddle-Node, Hopf, and Turing bifurcations’ boundary conditions are also examined for the system, and subsequently, Turing space is identified. The developed numerical scheme is applied to various multidimensional Gray–Scott systems to assess its ability to capture Turing pattern formation. Several test problems, such as stationary and non-stationary waves, pulse-splitting waves, self-replicating waves, and different types of spots, are chosen from the literature to demonstrate the different Turing patterns.