I have the following ODE for a function $y(x)$, over the real variable $x$, with $m>0$ an integer and $\epsilon$ a small parameter
$$ y'' + \frac{1}{x}y' + \frac{1}{x^{2}} y (\epsilon^{2} x^{2} -m^{2})=0, $$ with a requirement that the solution is finite at the origin ($x=0$).
Now, when $\epsilon=0$, the equation is reduced to the Cauchy-Euler equation, whose acceptable solution will be $c_{1}x^{m}$, where $c_{1}$ is a constant and we have rejected the other solution ($x^{-m}$) since it becomes unbounded at the origin. On the other hand, when $\epsilon\neq 0$, the equation becomes the parameterized Bessel equation, whose solution will be $c_{2} J_{m}(\epsilon x)$, where $c_{2}$ is a constant and $J_{m}$ is Bessel's function of the first kind and $m^{\textrm{th}}$ order. We have again rejected the second solution (second type Bessel function) as its magnitude shoots up to infinity at the origin.
Now, I wish to solve this as a perturbation problem in the small parameter $\epsilon$, by writing the solution as $$ y=y_{0}+\epsilon y_{1}+\epsilon^{2} y_{2}+\cdots.$$
When I do so by substitution into the original ODE, this gives me $y_{0}=c_{1}x^{m}$ and for $y_{1}$ the perturbed equation becomes an inhomogeneous Bessel equation as $$ y''_{1} + \frac{1}{x}y'_{1} + \frac{1}{x^{2}} y_{1} (\epsilon^{2} x^{2} -m^{2})+c_{1}\epsilon x^{m}=0, $$ whose solution is $y_{1}=y_{h}+y_{p}$, where $y_{h}$ denotes the associated homogeneous-equation solution and $y_{p}$ is a particular solution. We already know that $y_{h}=c_{2}J_{m}(\epsilon x)$, and for the particular solution, I tried to guess with trial polynomials and I found a suitable polynomial: $y_{p}=-c_{1} x^{m}/\epsilon$. So now we have $y_{1}=c_{2}J_{m}(\epsilon x)-c_{1} x^{m}/\epsilon$.
The problem now is when I want the full solution up to first order, I get $$ y\cong y_{0}+\epsilon y_{1}= c_{1}x^{m}+\epsilon[c_{2}J_{m}(\epsilon x)-c_{1} x^{m}/\epsilon] = \boxed{c_{2}\ \epsilon \ J_{m}(\epsilon x)}.$$
I am not sure if this makes sense because $y_{0}$ cancels out, and now the unperturbed problem gives 0 instead of $y_{0}$. I get the same result when I expand in powers of $\epsilon^{2}, \epsilon^{4}, \cdots$, instead of powers of $\epsilon, \epsilon^{2}, \cdots$ as above.
I thought about the fact that, for small argument, the Bessel function will become asymptotically $J_{m}(\epsilon x) \rightarrow \epsilon^{m} x^{m}2^{-m}/m!$, and therefore the result in the box above will be equivalent to $ y \rightarrow c_{2} \epsilon^{m+1} x^{m}2^{-m}/m!$, but again this gives zero for the unperturbed case, which doesn't seem to make sense!
What am I missing here?
Generally I find the habit of assuming that the powers of $\epsilon$ are in a particular form immediately is bad, because of how often problems arise where it doesn't happen. So generally I would prefer to write $y=y_0+\epsilon^\alpha y_1 + \dots$ and find $\alpha$ during the analysis.
Anyway, as you say, $y_0=c_1 x^m$; since everything is linear anyway we can choose a nontrivial IVP to solve, e.g. $y(1)=1$, which results in $c_1=1$. Now continue:
$$(x^m+\epsilon^\alpha y_1)''=m(m-1) x^{m-2}+\epsilon^\alpha y_1'' \\ \frac{1}{x}(x^m+\epsilon^\alpha y_1)'=m x^{m-2}+\frac{1}{x} \epsilon^\alpha y_1' \\ \frac{1}{x^2} (x^m+\epsilon^\alpha y_1)(\epsilon^2 x^2-m^2) = \frac{1}{x^2}(\epsilon^2 x^{m+2}-m^2 x^m+\epsilon^{2+\alpha} x^2 y_1-m^2 \epsilon^\alpha y_1).$$
The point is the term that you are looking to cancel out by including $y_1$ in the first place is $\epsilon^2 x^m$ in the last equation, so you should have $\alpha=2$ (not $1$). Then the equation for $y_1$ is
$$y_1''+\frac{1}{x} y_1' - \frac{m^2}{x^2} y_1 = -x^m$$
having discarded the $\epsilon^4$ term.