as follows, $$\begin{aligned}
\label{z3}
\frac{1}{\Omega}({\Phi}'_k+{\mathcal H}\Phi_k)=\frac{\bar{\phi}'}{2a}\left(1-\frac{k^2}{ (\bar{\phi}')^2  }\bar{f}_{\chi\chi}\right)\delta\sigma_k.\end{aligned}$$ Inserting $\Phi_k^{(0)}$ on it, one gets $\delta\sigma_k^{(0)}$. Then, to find the new iteration of $u_k^{(1)}$, one has to solve the equation $$\begin{aligned}
u_k''+k^2 u_k-\frac{\theta''}{\theta}u_k & = &-k^2\bar{\phi}'\sqrt{\Omega}\left( \frac{\bar{f}_{\chi\chi}\delta\sigma_k^{(0)}}{({\bar{\phi}')}^3 }  \right)',\end{aligned}$$ using the well-known [*method of variation of constants*]{} for second-order differential equations (see for example Chapts. $13-17$ of [@levinson]). When one has $u_k^{(1)}$ one calculates $\Phi_k^{(1)}$ and, inserting the result in (\[z3\]), one gets $\delta\sigma_k^{(1)}$. And, thus, the new iterations are iterativelly obtained and a fully fledged method is constructed.

Alternatively, one could also deal directly with Eq. (\[v\]), considering the right hand side as a perturbation. In Fourier space, at very early times this equation becomes $$\begin{aligned}
v''_k+k^2v_k-\frac{2}{\tau^2}v_k=0,\end{aligned}$$ whose solution satisfying the corresponding asymptotic is $$\begin{aligned}
\label{xxx}
v_k^{(0)}=\frac{e^{-ik\tau}}{\sqrt{2k}}\left(1-\frac{i}{k\tau}\right).\end{aligned}$$ When the pivot scale leaves the Hubble radius one can use the longwavelenght approximation (for the matter bounce scenario using the background (\[background\]) the function $\frac{z''}{z}$ is symmetric with respect to the expanding and contracting phase, being an increasing function in the whole contracting period. Thus, if the pivot scale leaves the Hubble radius at $t=-t_H$ it reenters at $t=t_H$) $$v''_k-\frac{z''}{z}v_k=0,$$ whose solution is [@haro13] $$\begin{aligned}
v^{(0)}_k=B_1(k)z+B_2(k)z\int_{-\infty}^{\tau}\frac{1}{z^2}d\bar\tau.\end{aligned}$$ At early times, if one deals 