Solve $y'' + \epsilon y'^2 + 1 = 0$ with initial conditions $y(0) = 0$ and $y'(0) = 1$ Let $\epsilon \ll 1$. 
I guess I'm trying to use perturbation method but I've been getting really weird numbers when I'm determining the initial conditions. Can someone perhaps help me with this?  
So what I have is the following:


*

*Set $y(t)=y_0(t)+\epsilon y_1(t)+\epsilon^2 y_2(t)+...$

*Plug in $y(t)$ into $y'$ and compare $\epsilon^{0}, \epsilon^{1}, \epsilon^{2}$:
$$(y_0''+\epsilon y_1''+\epsilon^{2}y_2''+...)+\epsilon(y_0'+\epsilon y_1'+\epsilon^{2}y_2'+...)^2+1=0$$
Next we collect all terms multiplied by $\epsilon$.  What I am doing here is basically finding all the coefficients of the $\epsilon$ and putting them together in a group: 
\begin{alignat}3
\epsilon^{0}&:& y_0''+1&=0& ~~ \text{where} ~~ y_0(0)&=0& ~~ \text{and} ~~ y_0'(0)&=1
\\
\epsilon^{1}&:& y_1''+y_0'^2&=0& ~~ \text{where} ~~ y_1(0)&=0& ~~ \text{and} ~~ y_1'(0)&=1
\\
\epsilon^{2}&:&~ y_2''+2y_0'y_1'+y_1'&=0& ~~ \text{where} ~~ y_2(0)&=0& ~~ \text{and} ~~ y_2'(0)&=1
\end{alignat}
I think I would treat them like regular ODEs but the numbers I'm getting are a bit weird.  Can anyone help me out?
 A: The equation to solve is, as inferred from the $\epsilon^1$ equation,
$$
y''+ϵy'^2+1=0.
$$
Apart from the fact that this is solvable using separation of variables or some Riccati approach, one can also treat this, as carried out, as a perturbation problem setting $$y=y_0+ϵy_1+ϵ^2y_2+...$$
Then inserting into the ODE one gets by comparing equal powers of $ϵ$
\begin{align}
y_0''+1&=0 &\implies y_0(x)&=x-\tfrac12x^2\\
y_1''+y_0'^2&=0&\implies y_1(x)&=-\tfrac1{12}(x^4-4x^3+6x^2)\\
y_2''+2y_0'y_1'&=0&\implies y_2''(x)&=\tfrac23(3x-3x^2+x^3)(1-x)=\tfrac23(3x-6x^2+4x^3-x^4)\\
&&y_2(x)&=\tfrac23(\tfrac12x^3-\tfrac12x^4+\tfrac15x^5-\tfrac1{30}x^6)
\end{align}
The initial condition are already accounted for in $y_0$, the other components have homogeneous initial conditions. In the quadratic terms of the equation for $ϵ^k$ the indices inside a term have to add to $k-1$, thus for $k=2$ there can only be the index combination $(0,1)$.

To test, compare with the Riccati approach that sets $y'=\frac{u'}{ϵu}$ with $u(0)=1$, $u'(0)=y'(0)ϵu(0)=ϵ$, so that in turn $y(x)=\frac1ϵ\ln|u(x)|+y(0)$, to get the linear equation $$u''+ϵu=0$$ which has the solution
$$
u(x)=\cos(\sqrtϵx)+\sqrtϵ\sin(\sqrtϵx)\\\implies y(x)=\tfrac1ϵ\ln(\cos(\sqrtϵx)+\sqrtϵ\sin(\sqrtϵx))
$$
Expanding into power series in $x$ gives then
\begin{align}
y(x)&=\tfrac1ϵ\ln(1+ϵx-\tfrac12ϵx^2-\tfrac16ϵ^2x^3+\tfrac1{24}ϵ^2x^4+...)
\\
&=(x-\tfrac12x^2-\tfrac16ϵx^3+\tfrac1{24}ϵx^4+\tfrac1{120}ϵ^2x^5-\tfrac1{720}ϵ^2x^6...)\\&~~~~-\tfrac12ϵ(x-\tfrac12x^2-\tfrac16ϵx^3+\tfrac1{24}ϵx^4+...)^2 + \tfrac13ϵ^2(x-\tfrac12x^2+...)^3+...
\\
&=x-\tfrac12x^2+ϵ(-\tfrac12x^2+\tfrac13x^3-\tfrac1{12}x^4)+ϵ^2(\tfrac13x^3-\tfrac13x^4+\tfrac2{15}x^5-\tfrac1{45}x^6)+...
\end{align}
which is as expected exactly the same as the perturbation series.
