How to find asymptotic expansions of all real roots of $x \tan(x)/\epsilon=1?$ 
Find expansions of all the real roots of  $$x\tan(x)=\epsilon?$$ (You have to consider the first root separately)

It is really bothering me. If I assume $x=x_0+x_1\epsilon +x_2\epsilon^2$ and do the taylor expansion of tan(x).
Then I end up with $$(x_0+x_1\epsilon +x_2\epsilon^2)[(x_0+x_1\epsilon +x_2\epsilon^2)+1/3 (x_0+x_1\epsilon +x_2\epsilon^2)^3+...]=\epsilon$$.
From the above equation $x_0=0$, then no matter what value $x_1$ is, the coefficient of $\epsilon$ on the LHS is always $0$. 
Could anyone kindly give me some hints? Thanks!

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.)
