Not yet reviewed by Pith; the record is open.

Perturbation theory is a powerful tool for studying large-scale structure formation in the universe and calculating observables such as the power spectrum or bispectrum. However, beyond linear order, typically this is done by assuming a simplification in the time-dependence of gravitational-coupling kernels between the matter and velocity fluctuations. Though the true dependencies are known for Lambda cold dark matter cosmologies, they are ignored due to the computational costs associated with considering them in full and, instead, are replaced by simpler dependencies valid for an Einstein--de-Sitter cosmology. Here we develop, implement and demonstrate the effectiveness of a new numerical method for finding the full dynamical evolution of these kernels to all perturbative orders based upon spectral methods using Chebyshev polynomials. This method is found to be orders of magnitude more efficient than direct numerical solvers while still producing highly accurate and reliable results. A code implementation of the Chebyshev spectral method is then presented and characterised. The code has been made publicly available alongside this paper. We expect our method to be of use for interpretation of upcoming galaxy clustering measurements.