Question about phase shift on multiple-scale analysis Consider the following ODE
$$y''(t) + y(t) + \epsilon y^2(t) y'(t) = 0$$
for $t>0$ with boundary condition $y(0)=1$ and $y'(0)=0$
I have found the leading order asymptotic expansion, that is
$$y(t)={2 \over \sqrt{\epsilon t +4}} \cos(t)$$
My question is why for $\epsilon$ small, like $0.1$, this asymptotic expansion matched pretty well with the numerical solution. However, if I increase $\epsilon$ to say $5$, they don't match at all. What I notice is that the amplitude of those graphs still match, just they are like out of phase. 
Here is the graph, orange one is asymptotic expansion and blue one is numerical solution.

 A: Generally, the formulation of a perturbation problem supposes that $ε$ is small relative to the other constants in the equation. Here these constants are $\sim 1$, so that $ε\ll 1$ is implied. This is not the case for $ε=5$

Your first order approximation also seems incomplete.
Setting $y(t)=r(t)·\cos(φ(t))$ and $y'(t)=-r(t)·\sin(φ(t))$ with, as per the zero-th order solution, $r=1+O(ε)$ and $φ=t+O(ε)$, one obtains
\begin{align}
r^2&=y^2+y'^2\\
rr'&=yy'+y'y''=-ϵy^2y'^2
&\implies 
\left(\frac1{r^2}\right)'&=-\frac{2r'}{r^3}=\frac{ε}4(1-\cos(4t))+O(ε^2)
\\ \\
\tan(φ)&=-\frac{y'}{y}\\
r^2φ'&=-y''y+y'^2=ϵy^3y'+r^2
&\implies
φ'&=1-εr^2·\cos^3(φ)·\sin(φ)-1\\
&&&=\left(t+\frac{ε}4\cos^4(t)\right)'+O(ε^2)
\end{align}
which integrates in the first order of $ε$ to
\begin{align}
\frac1{r^2}&=1+\frac{ε}4·\left(t-\frac14\sin(4t)\right)+O(ε^2)\\
φ&=t-\frac{ε}4\left(1-\cos^4(t)\right)+O(ε^2)
\end{align}
which combines to
$$
y(t)=\frac{1}{\sqrt{1+\frac{ε}4(t-\frac14\sin(4t))}}·\cos\left(t-\frac{ε}4\left(1-\cos^4(t)\right)\right)
$$
or in a simplistic average
$$
y(t)=\frac{1}{\sqrt{1+\frac{ε}4 t}}·\cos\left(t-\frac{5ε}{32}\right).
$$
This approximation still looks nice relative to the numerical solution for $ε=2$, the case $ε=5$ looks a little better in phase, but still shifts away for larger values of $t$.
