  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.
