Numerical solution to Mukhanov-Sasaki equation I am trying to figure out how to solve the Mukhanov-Sasaki equation to compute the power spectrum of an inflation potential that exhibits an ultra slow-roll phase that gives rise to an enhancement in the power spectrum, suitable for primordial black hole formation. In terms of the efold variable, the MS equation is
$$ \frac{d^{2}\zeta_{k}}{dN^{2}}+\left( 3-\epsilon_{1}+\epsilon_{2} \right)\frac{d\zeta_{k}}{dN}+\left( \frac{k}{aH} \right)^{2}\zeta_{k}=0 $$
with the Bunch-Davies vacuum as the initial condition
$$ u_{k}=\lim_{k>>aH}\frac{e^{-ik\tau}}{\sqrt{2k}} $$
I have computed the Hubble flow parameters but now I am unsure how to solve this ODE for all the modes of interest. I know that I need to consider them deep inside the horizon and evolve them to horizon exit where they freeze over, however, I am unsure how to proceed.
 A: You need to solve the MS equation separately for each $k$, initialized far enough in the past so that each mode is approximately in the Bunch-Davies limit.  I've found that choosing $N_i$ so that $k = \mathcal{O}(100)a(N_i)H(N_i)$ is sufficient for this.  Note that there is a different $N_i$ for each $k$.  Here, $a(N) = a_0 \exp(N_0-N)$, where you need to pick a value for $a_0$.  For example, you could set $a_0$ such that a particular scale, like the quadrupole, is at horizon crossing at $N_0$, i.e. $k_{\ell = 2} = a_0 H_0$.  
Initializing further back in time will help with accuracy, but there are diminishing returns.  Also, if you go back too far, you'll begin to see trans-Planckian modulations in your power spectra.
Then, you just evolve each mode forward in time until inflation ends.  The largest-scale modes will have frozen-out long before the end of inflation, but this way you can compute $P(k)$ by evaluating each individual mode $u_k$ on an equal-time slice.  Alternatively, you could evolve each mode until they are well-outside the horizon, say, when $k < aH/100$ or something.  The modes freeze-out quite soon after horizon exit so applying a generous cut-off like this should work well.  
