Your analysis is correct.
Let's look at the original equation
$${\epsilon ^2}{x^4} - \epsilon {x^3} - 2{x^2} + 2 = 0$$
and plot the Newton-Kruskal diagram:
There are only two possible placements of straight lines passing through two or more points with all remaining points "above the line"; those two lines correspond exactly to
- $\delta^2 \sim 2 \Rightarrow \delta \sim 1$, and
- $\epsilon^2\delta^4 \sim \epsilon \delta^3 \sim \delta^2 \Rightarrow \delta \sim \epsilon^{-1}.$
Before we continue, let's rescale the original equation by making the change of variables $x=\delta y$ to arrive at
$${\epsilon ^2}{\delta ^4}{y^4} - \epsilon {\delta ^3}{y^3} - 2{\delta ^2}{y^2} + 2 = 0.$$
From the first balance, $\delta^2 \sim 2 \Rightarrow \delta \sim 1$, we immediately recover the original equation:
$$\tag{1} {\epsilon ^2}{y^4} - \epsilon {y^3} - 2{y^2} + 2 = 0.$$
Letting $\epsilon \rightarrow 0$, we get
$$2{y^2} + 2 = 0 \Rightarrow y = \pm 1$$
Let´s check if there are roots close to $y=-1$ and $y=1$.
Plugging an asymptotical expansion of the form $y \sim \pm 1 + {a_1}\epsilon + {a_2}{\epsilon ^2} + {a_3}{\epsilon ^3} + {a_4}{\epsilon ^4} + \ldots $ into (1), and then matching coefficients, we arrive at
$${x_a} \sim {y_a} \sim - 1 - \frac{\epsilon }{4} - \frac{{13}}{{32}}{\epsilon ^2} + O\left( {{\epsilon ^3}} \right),$$
$${x_b} \sim {y_b} \sim 1 - \frac{\epsilon }{4} + \frac{{13}}{{32}}{\epsilon ^2} + O\left( {{\epsilon ^3}} \right).$$
Regarding the second possibility of balancing, $\epsilon^2\delta^4 \sim \epsilon \delta^3 \sim \delta^2 \Rightarrow \delta \sim \epsilon^{-1}$, we get the following equation
$$\tag{2}{y^4} - {y^3} - 2{y^2} + 2{\epsilon ^2} = 0$$
Letting $\epsilon \rightarrow 0$, we get
$${y^4} - {y^3} - 2{y^2} = 0 \Rightarrow {y^2}\left( {{y^2} - y - 2} \right) = 0 \Rightarrow y = 0 \vee y = - 1 \vee y = 2.$$
Let´s check if there are roots close to $y=-1$ and $y=2$ ($y=0$ is of no interest, as we just recover one of the roots previously found). Plugging an asymptotical expansion of the form $y \sim {a_0} + {a_1}\epsilon + {a_2}{\epsilon ^2} + {a_3}{\epsilon ^3} + {a_4}{\epsilon ^4} + \ldots $ (with ${a_0} = - 1 \vee {a_0} = 2$) into (2), and then matching coefficients, we arrive at
$${x_c} \sim {\epsilon ^{ - 1}}{y_c} \sim - {\epsilon ^{ - 1}} + \frac{2}{3}\epsilon + O\left( {{\epsilon ^3}} \right),$$
$${x_d} \sim {\epsilon ^{ - 1}}{y_d} \sim 2{\epsilon ^{ - 1}} - \frac{1}{6}\epsilon + O\left( {{\epsilon ^3}} \right).$$
Thus, we've determined all the four roots of the original equation.
-EDIT-
Here's the Mathematica code I used to do my calculations, more specifically, the last one.
n = 4;
a[0] = 2;
x = a[0] + Sum[a[i] \[Epsilon]^i, {i, 1, n}] + O[\[Epsilon]]^(n + 1);
x^4 - x^3 - 2 x^2 + 2 \[Epsilon]^2 == 0;
LogicalExpand[%]
Solve[%]
Best Answer
Your approach is working just fine.
At $O(1)$ you have
$$ x_0^3-x_0^2=-0$$
with solutions $x_0=0,0,1$. Then at $O(\epsilon)$ your equation depends on the value of $x_0$. If $x_0=1$, you get
$$ 3x_1-2x_1-1=0\Rightarrow x_1=1. $$
If $x_0=0$ then the next equation you get is at $O(\epsilon^2)$, and
$$ -x_1^2+1=0\Rightarrow x_1=\pm1. $$
So it depends a bit on what you want for your result, because $x_0=1+\epsilon$ satisfies the equation to within $O(\epsilon)$, but $x_0=\pm\epsilon$ satisfy the equation to within $O(\epsilon^2)$ (notice that $x=0$ satisfies the original equation to within $O(\epsilon^2)$), but all three approximations at within $O(\epsilon)$ of the real root.