I started with the following 2nd order PDE $$ -\frac{x}{x-1}\partial_y^2\phi+\frac{1}{x^2} \partial_x (x(x-1)\partial_x \phi) = \mu^2 \phi$$ which is clearly separable. I want a solution that asymptotes to a $y$-dependent constant at $x$-infinity. I therefore tried the form $\phi = X(x)Y(y)$, with $X(x\rightarrow \infty) \rightarrow1$. Looking at this limit above we see that it can only hold if $\partial_y^2 Y = -\mu^2 Y$. Subbing this back in, I end up with an ODE for $X$: $$ \frac{1}{x^2} \partial_x (x(x-1)\partial_x X) = -\frac{\mu^2}{x-1} X$$ I've tried solving this numerically (obviously staying away from $x=1$) on Mathematica. For any initial conditions I give at some finite $x$, the solution seems to always die off oscillating around zero, as in plot.
So I'm wondering why it seem like the boundary condition at infinity $X\rightarrow $ constant seems allowed. Have I fooled myself?
If I expand the ODE for large $x$ to 'first order' I have $$\left( 1- \frac{1}{x}\right)\partial_x^2 X + \frac{2}{x} \partial_x X = -\frac{\mu^2}{x} X, \quad \quad x \rightarrow \infty$$ Mathematica gives exact solutions for this in terms of Bessel functions and these indeed match the asymptotic behaviour in the plot above. However if I stopped at 'zeroth order' I would just have $\partial_x^2 X = 0$ and a self-consistent solution $X=$ constant (a term linear in $x$ would not be self-consistent with ignoring the first derivative term).
Where is my reasoning wrong? Is there no solution with the boundary condition I desire? (If there is a way of enforcing a boundary condition at infinity numerically I am unaware of it)
First, let me very briefly restate the argument you have presented. "If we assume $\lim_{x \rightarrow \infty} X(x) = 1$, then we discover that $\lim_{x \rightarrow \infty} X(x) = 0$, which is a contradiction."
With a little more work, if we assume $\lim_{x \rightarrow \infty} X(x) = k$, for some constant $k$, then we have (starting in the middle and evaluating the two limits as we go out) $$ k Y''(y) = \lim_{x \rightarrow \infty} \frac{k x Y''(y)}{x-1} = \lim_{x \rightarrow \infty} k \mu^2 Y(y) = k \mu^2 Y(y) \text{,} $$ obtaining your $Y'' = \mu^2 Y$ or $k = 0$. (If $k \neq 0$, the division to cancel the $k$s is valid. If $k = 0$, the division is invalid. We have to handle each case separately.) You have shown that if we assume $k \neq 0$, we get the contradiction $k = 0$. So now we consider $k= 0$, but there is no information in the equation $0 \cdot Y''(y) = 0 \cdot \mu^2 Y(y)$.
What we conclude from the argument you have given is : if $\lim_{x \rightarrow \infty} X(x)$ exists, that limit is zero; however this gives no constraint on the differential equation in the limit.
Generically, a way to investigate the behaviour around $x \rightarrow \infty$ is to make the change of variable $x \mapsto 1/u$. Then $\partial_x = \frac{\mathrm{d}x}{\mathrm{d}u} \partial_u = -u^{-2} \partial_u$. However, we can do this most expediently at the point of declaring your separation of variables: $\phi = U(1/x)Y(y)$ and then replace every $x$ with a $1/u$. Then consider $u \rightarrow 0^+$. After substitution and a little manipulation, (where primes are derivatives with respect to slots not variables) $$ \frac{ (u-1) (\mu^2 U(u) + u^4( U'(u) + (u-1)U''(u))}{U(u)} = \frac{Y''(y)}{Y(y)} \text{.} $$ Taking $u \rightarrow 0^+$, we have $$ \lim_{u \rightarrow 0^+} -\left( \mu^2 + u^4\frac{U'(u) + (u-1)U''(u)}{U(u)} \right) = \frac{Y''(y)}{Y(y)} \text{.} $$ And now we have to try to extract information about the order of zero of $U$ at $u = 0$. Which is hard. (As an example of why this is not easy, take $U(u) = \mathrm{e}^{-1/u^2}$, so that the limit of the left hand side is $\infty$, not $-\mu^2$. The point is not that $\mathrm{e}^{-1/u^2}$ is a solution, it is that determining the algebraic multiplicity of zeroes and poles in functions given by differential equations need not be easy.)
Alternatively, let's use what we suspect about $Y(y)$, specifically, $\frac{Y''(y)}{Y(y)} = - \mu^2$, so that the ansatz is $Y(y) = \mathrm{e}^{\mathrm{i} \mu y}$ and we take $\phi(x,y) = X(x) \mathrm{e}^{\mathrm{i} \mu y}$, giving $$ \mathrm{e}^{\mathrm{i} \mu y} \left( \frac{x \mu^2 X(x)}{x-1} + \frac{(2x-1) X'(x) + (x-1)x X''(x))}{x} \right) = 0 \text{,} $$ which, since the exponential is never zero, forces $$\frac{x \mu^2 X(x)}{x-1} + \frac{(2x-1) X'(x) + (x-1)x X''(x))}{x} = 0 \text{.} $$ Trying $x \rightarrow \infty$, we find $$ \mu^2 X(x \rightarrow \infty) + 2 X'(x \rightarrow \infty) + \lim_{x \rightarrow \infty} (x-1)X''(x) = 0 \text{.} $$ We can solve this and discover that asymptotically, $X(x)$ behaves like $$ \frac{c_1 \sqrt{\mu^2 (1-x)} I_1(2 \sqrt{\mu^2(1-x)})}{1-x} - \frac{c_2 \sqrt{\mu^2 (1-x)} K_1(2 \sqrt{\mu^2(1-x)})}{1-x} \text{.} $$ From this, we see that if $c_1 = 0$, $|X| \rightarrow 0$ and if $c_1 \neq 0$, $|X| \rightarrow \infty$. Thus, (assuming our ansatz) if $X$ approaches some constant as $x \rightarrow \infty$, that constant is zero, and generally $X$ need not be bounded as $x \rightarrow \infty$.
Maybe some other analysis lets us see that $c_1 = 0$, but I don't know what that analysis is.