Rapid numerical solutions for the Mukhanov-Sasaki equation

We develop a novel technique for numerically computing the primordial power spectra of comoving curvature perturbations. By finding suitable analytic approximations for different regions of the mode equations and stitching them together, we reduce the solution of a differential equation to repeated matrix multiplication. This results in a wavenumber-dependent increase in speed which is orders of magnitude faster than traditional approaches at intermediate and large wavenumbers. We demonstrate the method's efficacy on the challenging case of a stepped quadratic potential with kinetic dominance initial conditions.