Recently, I have received questions regarding the small $\epsilon$ expansion for the largest real root of
$$\tag{1} f(y)=\epsilon y (y^2-1)^2-(2y^2-1) $$
I am posting this self-answered question as a reference for the askers. The direct perturbation ansatz $y=\sum\limits_n\epsilon^ny_n$ fails for the largest root because as $\epsilon \to 0$ this root goes to infinity, while at $\epsilon=0$ there are only two roots. This is a singular perturbation problem that arises in the study of the equilibrium position of a charge near a conducting sphere (see eg p61 in Jackson).
The first step in a singular perturbation problem is to rescale $y$ such that the unknown is of order unity as $\epsilon\to0$. Let
$$\tag{2} y(\epsilon)=x(\epsilon)\epsilon^\alpha $$
then substituting (2) into (1) yields the relation
$$\tag{3} x^3\epsilon^{5\alpha+1}\sim2\epsilon^{2\alpha} \qquad ,\qquad \epsilon\to0 $$
so we choose $\alpha=-1/3$ and immediately read off the leading order behavior $x\sim2^{1/3}$ as $\epsilon\to0$. We need to find the approximate largest root of
$$\tag{4} F(x)=\epsilon^{-2/3}\left(x^5-2x^2\right)+(1-2x^3)-\epsilon^{2/3}x $$
Using (4), one readily verifies that $x(\epsilon)$ does not have an expansion in integral powers of $\epsilon$. In fact, the consistent series is
$$\tag{5} x(\epsilon)=\sum\limits_{n=0} \epsilon^{2n/3}x_n $$
It is then straightforward to determine the coefficients $x_n$ recursively. The first few terms in the expansion of $y$ are
$$\tag{6} y(\epsilon)\sim \left(\frac{\epsilon}{2}\right)^{-1/3} +\frac{1}{2}\left(\frac{\epsilon}{2}\right)^{1/3} -\frac{1}{12}\left(\frac{\epsilon}{2}\right)^1+\frac{1}{36}\left(\frac{\epsilon}{2}\right)^{7/3}-\frac{1}{32}\left(\frac{\epsilon}{2}\right)^3\quad,\quad \epsilon\to0 $$
This is nice, but we can do better. Taking the Padé approximants $P_M^N$ of (6) yields approximations that appear valid all the way to large $\epsilon$. The first diagonal Padé approximant is
$$\tag{7} y_1^1(\epsilon)=\left(\frac{\epsilon}{2}\right)^{-1/3}\left[\frac{6+4(\epsilon/2)^{2/3}}{6+(\epsilon/2)^{2/3}}\right] $$
In the plot below we compare the numerically found root (black line) to the partial sum and to $y_1^1$ for 'small' values of $\epsilon<5$ (dashed lines). The dotted lines are the percentage error of each approximation.
The same procedure applied to the large $\epsilon$ regime yields
$$\tag{8} y(\epsilon)\sim 1 + \frac{1}{2\sqrt{2}}\left(\frac{\epsilon}{2}\right)^{-1/2}+\frac{1}{8}\left(\frac{\epsilon}{2}\right)^{-1}\quad,\quad \epsilon\to\infty $$
Using (6) and (8) we can form a two-point Padé approximant using the points $\epsilon=0$ and $\epsilon=\infty$. The simplest such approximation is
$$\tag{9} y^{\text{approx}}=\frac{\beta-6}{\beta^2+2\beta+6}+1+\frac{1}{\beta} $$
where $\beta=(\epsilon/2)^{1/3}$.