Perturbation method to solve $\epsilon x^3=x+\epsilon$ I was given the following exercise:

Find two-term approximation to the real roots of the equation
$\epsilon x^3=x+\epsilon \tag*{}$

My attempt:
Since this is a singular equation, we introduce a change of variable $x=\delta y$ with $y\simeq 1$. Substitute this to the equation and we get
$\epsilon \delta^3 y^3 = \delta y +\epsilon \tag*{}$
To find the appropriate form of $\delta$, we use the method of dominant balance. We take the balance $\epsilon \delta^3 y^3\simeq \delta y$ and we have $\delta \simeq \frac{1}{\sqrt{\epsilon}}$. Substitute back and we get
$\frac{y^3}{\sqrt{\epsilon} } = \frac{y}{\sqrt{\epsilon}} +\epsilon\tag*{}$
or
$y^3 =y +\epsilon^{3/2}\tag*{}$
Expand $y$ with $y=y_0 +y_1 \epsilon+\cdots $. Then we have
$y_0^3+3y_0^2y_1\epsilon +... = y_0+y_1\epsilon +\epsilon^{3/2}\tag*{}$
If we compare the coefficient of $\epsilon^{3/2}$ term, we have $0=1$.
(I am stuck here)

Since the expansion $y=y_0 +y_1 \epsilon+...$ gives contradiction, I think we need to use another expansion, but I am not sure which kind of expansion I should use.
 A: Let's start from $y^3=y+\epsilon^{3/2}$ from the question.
Let $y = y_0 + y_1\epsilon^{\gamma}+\cdots$ where $\gamma$ is a constant.
$y_0$ is the solution to $y^3=y$. That is, $y_0=0,1,-1$
We use the method of dominant balance to find the appropriate $\gamma$. By substituting the expansion to the equation, we get the following equation.
$y_0^3 + 3y_0^2y_1\epsilon^{\gamma}+\cdots = y_0 + y_1 \epsilon^{\gamma}+\cdots + \epsilon^{3/2} \tag*{}$
Since $y_0^3=y_0$,
$3y_0^2y_1\epsilon^{\gamma}+\cdots =y_1 \epsilon^{\gamma}+\cdots + \epsilon^{3/2} \tag*{}$
We consider $\epsilon^{\gamma}\simeq \epsilon^{3/2}$ and every other term is negligible. Then we have $\gamma \simeq 3/2$. We can easily see that this value is consistent with the balance. Therefore set $\gamma = 3/2$.
Now consider $O(\epsilon^{3/2})$ terms. By comparing the coefficient, we have
$3y_0^2y_1=y_1+1\tag*{}$
By solving this, we have
$\displaystyle y_1 = \frac{1}{3y_0^2-1}\tag*{}$
For $y_0=0,1,-1$, the value is $y_1 =-1, 1/2, 1/2$, respectively. Therefore, we have three asymptotic solution for $y$:
$\begin{align} y &= -\epsilon^{3/2}+\cdots \\ y &= 1 +\frac{\epsilon^{3/2}}{2} +\cdots \\ y=&= -1 +\frac{\epsilon^{3/2}}{2} + \cdots \\ \end{align}\tag*{}$
By using $x=y/\sqrt{\epsilon}$, we have the following solutions:
$\begin{align} x &= 0-\epsilon +\cdots \\ x &= \epsilon^{-1/2} +\frac{\epsilon}{2} +\cdots \\ x &= -\epsilon^{-1/2} +\frac{\epsilon}{2} + \cdots \\ \end{align}\tag*{}$
Observe that one of the solutions matches the expansion mentioned in the comment.
