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.

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.