  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.
