Iterative method for finding asymptotic expansion of solution Question:
For the cubic equation
$$x^3+x-\varepsilon = 0 \qquad \varepsilon \rightarrow 0$$
find the first three terms of the asymptotic expansion of the solutions.

Attempt:
I observed that a solution is $x=0$ when $\varepsilon = 0$, so I picked the starting point $x_0=0$.
I then rearranged the equation as $x = (\varepsilon - x)^{1/3}$, to give the recurrence relation
$$x_{n+1} = (\varepsilon - x_n)^{1/3}$$
Plugging it in, we get
\begin{align}
x_0 & = 0 \\
x_1 & = \varepsilon^{1/3} \\
x_2 & = (\varepsilon - \varepsilon^{1/3})^{1/3} \sim \mathcal O(\varepsilon^{1/9})
\end{align}
and I am confused, because apparently future iterations are supposed to match previous iterations at the leading orders.
But in this case, I had $x_1 = \mathcal O(\varepsilon^{1/3})$ and $x_2 = \mathcal O(\varepsilon^{1/9})$.
What have I done wrong? Any hints are much appreciated. Please help T_T
 A: The behaviour you are expecting will be evident if you switch to a parameter going to infinity.  So let $\nu = 1/\varepsilon$.  Then $\nu \rightarrow \infty$ as $\varepsilon \rightarrow 0$.
For instance, at the step you are calling out in your computation, 
$$  (\nu^{-1} - \nu^{-1/3})^{1/3} = \nu^{-1/3} - \frac{1}{3} \nu^{1/3} - \frac{1}{9} \nu - \cdots  \text{.}  $$
(Some people also just write as if they are expanding in $\frac{1}{\varepsilon}$.)
A: By the Lagrange inversion theorem, in a neighbourhood od the origin the inverse function of $f(x)=x+x^3$ can be written as
$$ f^{-1}(x) = \sum_{n\geq 1}\frac{x^{n}}{n}\cdot [z^{n-1}]\left(\frac{1}{1+z^2}\right)^n=\sum_{k\geq 0}\frac{x^{2k+1}}{2k+1}\cdot [z^{2k}]\left(\frac{1}{1+z^2}\right)^{2k+1} $$
and by stars and bars
$$ [z^{2k}]\left(\frac{1}{1+z^2}\right)^{2k+1}=[z^{k}]\left(\frac{1}{1+z}\right)^{2k+1}=(-1)^k[z^k]\left(\frac{1}{1-z}\right)^{2k+1}=(-1)^k\binom{3k}{k} $$
so
$$ f^{-1}(x) = \sum_{k\geq 0}\frac{(-1)^k x^{2k+1}}{2k+1}\binom{3k}{k}$$
and
$$\boxed{ f^{-1}(\varepsilon) = \varepsilon - \varepsilon^3 + 3 \varepsilon^5 + O(\varepsilon^7).}$$
See also Bring radical.
A: Let $x=\epsilon^\alpha x_0+\epsilon^\beta x_1+\epsilon^\gamma x_2+\ldots$, with $\alpha<\beta<\gamma<\ldots$. Then
$$ \epsilon^{3\alpha}x_0^3+3\epsilon^{2\alpha+\beta}x_0^2x_1+\epsilon^{\alpha}x_0+\epsilon^\beta x_1+O(\epsilon^{\alpha+2\beta},\epsilon^\gamma)=\epsilon. $$
So either $\alpha=0$ or $\alpha=1$ are possible balances. $\alpha=0$ gives $x_0^3+x_0=0$, so $x_0=0$, and $\alpha=1$ gives $x_0=1$. So let $\alpha=1$, and then the equation will be
$$\epsilon^3+3\epsilon^{2+\beta}x_1+\epsilon^\beta x_1+O(\epsilon^{1+2\beta},\epsilon^\gamma)=0 $$
which gives $\beta=3$, and $x_1=-1$.
This can be continued to generate terms in the expansion $x\approx \epsilon-\epsilon^3+3\epsilon^5-12\epsilon^7+O(\epsilon^9)$.
To do this iteratively, use the other option,
$$x_{n+1}=\epsilon-x_n^3$$
which, starting with $x_0=0$, gives
$$ x_1=\epsilon, \quad x_2=\epsilon-\epsilon^3,\quad x_3=\epsilon-\epsilon^3+3\epsilon^5-3\epsilon^7+\epsilon^9, \ldots$$
and so on. This way you're increasing the power on $\epsilon$ each time you iterate, rather than reducing it.
