Dense colloidal suspensions exhibit pronounced viscoelasticity, shear‑thinning and yield‑stress behaviour arising from microscopic interactions and memory effects. Mode‑coupling theory (MCT) provides constitutive equations derived from first principles in statistical mechanics, resulting in an integral formulation for the definition of the extra stress tensor, along with integro‑differential equations for microstructural correlations. The nonlinear coupling between these fields poses severe numerical challenges typical of high‑Weissenberg‑number regimes. We formulate the problem within an operator‑splitting framework that decouples kinematic transport from memory evolution. Spatial discretisation is performed using discontinuous Galerkin methods, enabling local conservation and high parallel efficiency. A logarithmically spaced grid in the material's "age"' variable captures long‑time relaxation dynamics with efficiently allocated memory. The resulting solver is implemented in FEniCSx and enables predictive simulations of MCT‑based viscoelastic flows, including distinctive phenomena such as flow‑induced residual stresses. This work establishes a computational foundation for new benchmark problems in microscopic rheological modelling.