Your method is valid for the roots away from $x=0$, but you missed a bunch of roots of $\tan$ that lead to non-trivial solutions. If $\epsilon=0$, then the solutions are $x=k\pi$ for any integer $k$.
The case $x=0$ is special, because $x\tan(x)$ has slope 0 here, and hence a single solution for $\epsilon=0$, but as soon as $\epsilon$ is not zero, that root splits into two roots, one in $(-\pi/2,0)$ and the other in $(0,\pi/2)$, so it's a singular perturbation problem. I'll do this case after the $x\neq0$ case.
For the $x\neq 0$ case, you have $x=k\pi+\epsilon x_1+\epsilon^2x_2+\cdots$, and you know that $\tan(x)$ will be small because $\tan(k\pi)=0$ and $\tan$ is continuous. So, expand $\tan(k\pi+\epsilon x_1+\epsilon^2x_2+\cdots)$ around $k\pi$:
$$ \tan(k\pi+\epsilon x_1+\epsilon^2x_2+\cdots)=\tan(k\pi)+\frac{\epsilon x_1+\epsilon^2x_2+\cdots}{\cos^2(x_0)}-(\epsilon x_1+\epsilon^2x_2+\cdots)^2\frac{2\tan(x_0)}{\cos^2(x_0)}+O(\epsilon^3) $$
where I've gone up to $O(\epsilon^2)$. Putting in $x_0=k\pi$ and truncating $O(\epsilon^3)$ terms gives
$$ \tan(k\pi+\epsilon x_1+\epsilon^2x_2+\cdots)=\epsilon x_1+\epsilon^2x_2+O(\epsilon^3) $$
and now we are looking to solve
$$ (k\pi+\epsilon x_1+\epsilon^2x_2)(\epsilon x_1+\epsilon^2x_2)=\epsilon $$
the $O(\epsilon)$ equation is
$$k\pi x_1=1\Rightarrow x_1=\frac{1}{k\pi}. $$
Then the $O(\epsilon^2)$ equation is
$$x_1^2+k\pi x_2=0\Rightarrow x_2=-\frac{x_1^2}{k\pi}=-\frac{1}{(k\pi)^3}. $$
So for the roots away from $x=0$ you get $$x=k\pi+\frac{\epsilon}{k\pi}-\frac{\epsilon^2}{(k\pi)^3}+O(\epsilon^3).$$
For the roots near 0, we already determined it was a singular perturbation problem, so the standard expansion $x=x_0+\epsilon x_1+\epsilon^2x_2+\cdots$ wont work. Instead, let's search for an expansion of the for $x=\epsilon^\alpha(x_0+\epsilon x_1+\epsilon^2x_2+\cdots)$. Use the Taylor series expansion of $tan$ to give
$$ x\tan(x)=\epsilon^\alpha(x_0+\epsilon x_1+\epsilon^2x_2+\cdots)\left(\epsilon^\alpha(x_0+\epsilon x_1+\epsilon^2x_2+\cdots)+\frac{1}{3}\epsilon^{3\alpha}(x_0+\epsilon x_1+\epsilon^2x_2+\cdots)^3+O(\epsilon^{5\alpha})\right) $$
which can be simplified to
$$ x\tan(x)=\epsilon^{2\alpha}\left((x_0+\epsilon x_1+\epsilon^2x_2+\cdots)^2+\frac{1}{3}\epsilon^{2\alpha}(x_0+\epsilon x_1+\epsilon^2x_2+\cdots)^4+O(\epsilon^{4\alpha})\right) $$
Now since $x\tan(x)=\epsilon$, we write
$$ \left((x_0+\epsilon x_1+\epsilon^2x_2+\cdots)^2+\frac{1}{3}\epsilon^{2\alpha}(x_0+\epsilon x_1+\epsilon^2x_2+\cdots)^4+O(\epsilon^{4\alpha})\right)=\epsilon^{1-2\alpha} $$
now the largest terms of the left hand side are $O(1)$, so to ensure the leading order solution is non-trivial we need the left hand side to be $O(1)$, and so $1-2\alpha=0\Rightarrow\alpha=1/2$.
Now our equation to solve is
$$ \left((x_0+\epsilon x_1+\epsilon^2x_2+\cdots)^2+\frac{1}{3}\epsilon(x_0+\epsilon x_1+\epsilon^2x_2+\cdots)^4+O(\epsilon^{2})\right)=1 .$$
The $O(\epsilon^0)$ equation is
$$x_0^2=1\Rightarrow x_0=\pm1.$$
The $O(\epsilon)$ equation is $$2x_0x_1+\frac{1}{3}x_0^4\Rightarrow x_1=\pm\frac{1}{6}.$$
So the first two terms in solution for the roots near zero are
$$ x=\epsilon^{1/2}-\frac{1}{6}\epsilon^{3/2},\quad\text{and}\quad x=-\epsilon^{1/2}+\frac{1}{6}\epsilon^{3/2}.$$
So we have
$$x=k\pi+\frac{\epsilon}{k\pi}-\frac{\epsilon^2}{(k\pi)^3}+O(\epsilon^3),\quad k\in\mathbb Z,\quad k\neq0$$
and
$$ x=\epsilon^{1/2}-\frac{1}{6}\epsilon^{3/2}+O(\epsilon^{5/2}),\quad\text{and}\quad x=-\epsilon^{1/2}+\frac{1}{6}\epsilon^{3/2}+O(\epsilon^{5/2}).$$
(Barring any arithmetic errors I have probably made! Although what I've posted seems to be alright looking at some plots.)
$$f({\bf x})=\left.\sum_{n=0}^{\infty} \frac 1 {n!}\left(\sum_{i=1}^{m}(x_i-x_{i_0})\frac {\partial}{\partial x_i} \right)^n f({\bf x})\right|_{{\bf x}={\bf x}_{0}} $$
where $m$ is the number of independent variables.
In your case (around $(x,y)=(0,0)$), it simplifies to:
$$f(x,y)=\sum_{n=0}^{\infty} \frac 1 {n!}\left(x\frac {\partial}{\partial x}+y\frac {\partial}{\partial y} \right)^n f(x,y)|_{(x,y)=(0,0)} $$
$$=f(0,0)+(xf_x(0,0)+yf_y(0,0))+\frac 1 2(x^2f_{xx}(0,0)+2xyf_{xy}(0,0)+y^2f_{yy}(0,0)) \dots$$
Best Answer
Here's a hint from my instructor. Rearranging the given equation as follows \begin{align*} x^2 + \varepsilon & = (e^x - 1)^2 \\ \varepsilon & = (e^x - 1)^2 - x^2 \\ & = (e^x - 1 - x)(e^x - 1 + x) \end{align*} You can then Taylor expand $e^x$ around 0, substitute the given asymptotic expansion and try to balance term. Note that for $\varepsilon=0$, $x=0$ is the solution to the equation $1 + \sqrt{x^2} = e^x$, so you may choose $x_0 = 0$.