Summary: We introduce a new approach to the numerical approximation of the time-fractional Navier-Stokes-Fokker-Planck (TFNSFP) system, which involves the coupling of the incompressible Navier-Stokes equations to a time-fractional Fokker-Planck equation. The model arises in the context of dilute polymeric fluids, and it enhances the standard integer-derivative version of the model by including memory effects. The key challenge associated with the numerical solution of the TFNSFP system is that, in addition to it being nonlocal in time, it is, even in its simplest form, posed on a spatial domain that is the Cartesian product of two \(d\)-dimensional domains, for \(d\in\{2, 3\}\), the \(d\)-dimensional flow domain and the \(d\)-dimensional configuration space domain. We combine a kernel compression technique based on rational approximation of the integral kernel of the time-fractional derivative with a space-splitting method. By doing so, we transform the time-fractional partial differential equation into a fixed number of integer-order in time partial differential equations. The Fokker-Planck equation posed on the \(2d\)-dimensional domain is decoupled into two \(d\)-dimensional problems, and a standard combination of the Hermite spectral method on the configuration space domain and a finite element method on the flow domain is then applied. Finally, we combine our numerical scheme for the time-fractional Fokker-Planck equation (TFFPE) with a standard Navier-Stokes solver. We propose an efficient implementation strategy based on the algebraic structure of the discretized time-fractional Fokker-Planck equation, which significantly reduces the computational cost.