Radial Schrodinger equation with inverse power law potential
What they show in the paper is that, for $\beta>2$, there are no solutions with the given asymptotic form as $r\to 0$ unless we assume $\alpha>0$. I think you can go further and show that there are no nonsingular solutions for $\beta>2,\alpha<0$, but I'm not sure.
What does this mean physically? Well, when we have a potential with a singularity, usually we think of it as just an approximation that breaks down sufficiently close to the singularity. For instance, we model the hydrogen atom with a potential $V\propto r^{-1}$, but really the potential near $r=0$ doesn't go to infinity, due to the nonzero size of the nucleus. We get away with using the singular potential because the "bad" behavior at $r=0$ doesn't qualitatively change the solutions. (And of course people do correct for nonzero-nuclear-size effects in atomic physics.)
If it's true that the ground state of the Schrodinger equation is singular for potentials of the given form ($\alpha<0,\beta> 2$), what that means is that this procedure doesn't work. To be specific, suppose that you solved the Schrodinger equation for a potential that looks like the given one down to some "cutoff" $r_0$, and is constant for smaller $r$. What you'd find is that the solution doesn't tend to some limit as $r_0\to 0$ -- the solution depends qualitatively on the size of that cutoff, no matter how small it is.
To answer your last question, for any repulsive potential ($\alpha,\beta>0$), you expect to find only continuum (unbound) states. Those states have $E>0$, and all positive values of $E$ are allowed. So if you try to solve numerically for the ground state, I'm not surprised that you seem to get zero.
Addition: After the discussion in the comments, it occurs to me that we can see why the case $\beta=2$ behaves the way it does. The Schrodinger equation in that case is $$ -{1\over 2}\nabla^2\psi + {\alpha\over r^2}\psi=E\psi. $$ Suppose that you'd found a bound-state solution $\psi_0$ corresponding to some energy $E_0<0$. Define a new solution by simply scaling the radial coordinate: $$ \psi_1(r)=\psi_0(cr) $$ for any $c>0$. Then $\psi_1$ is also a solution to the Schrodinger equation, with energy $E_1=c^2E_0$. In particular, for $c>1$ this corresponds to squeezing the wavefunction into a smaller space and making the energy more negative. If you try to find the lowest-energy solution, you'll end up with the $c\to\infty$ case -- an infinitely concentrated wavefunction, with energy $-\infty$.
If you try the procedure I suggest in my last comment (cutting off the singularity in the potential at some $r_0$ and then varying $r_0$), something similar occurs. The ground state solution for all positive $r_0$'s look the same, with radial coordinates scaled by the value of $r_0$, and the ground state energy goes like $r_0^2$. As $r_0\to 0$, the ground-state energy approaches $-\infty$, and the wavefunction becomes infinitely concentrated at $r=0$.
This only works for the case $\beta=2$, because for this value of $\beta$ both the kinetic and potential terms on the left side of the Schrodinger equation scale in the same way when you rescale your coordinates (i.e., both go like $c^2$). Another way to put it: only in the case $\beta=2$ is the constant $\alpha$ dimensionless. For any other $\beta$, the value of $\alpha$ determines a length scale (so that you can't just rescale one solution to get a new one), but when $\beta=2$ the problem is scale-invariant.
This is addressing NGY' questions in the comment (and following up on my own comment). Suppose $\phi(r) \sim r^n$ for $r \sim 0$.
For $\beta < 2$ the derivate term beats the potential term and we get that around the origin the solution must (up to higher order terms) behave as $$n(n-1)r^{n-2} \sim 0$$ which happens for $n > 2$ and $n=1$. Note that if the solution were exactly $n=1$ (higher order terms included) then the derivative term would vanish identically and the equation would not be satisfied for $\beta \leq 1$. As a sidenote recall that there are also bound solutions that diverge for $r \to 0$ but they can still be normalizable, since normalization is given by an integral over the whole space and contribution from the sphere $S^{d-1}$ can beat the divergence from $\left | \phi(r)^2 \right |$.
For $\beta = 2$ we see that $n(n-1) r^{n-2} \sim \alpha 2 r^{n-2}$ implying $n=-1$ for $\alpha = 1$ and $n=2$ for $\alpha = -1$.
For $\beta > 2$ the potential term beats the derivative and we get $r^{n-\beta} = 0$. Which is satisfiable for any $n > \beta$ but then the solution will have to vanish identically at all orders in order to satisfy the equation (since the derivative, the potential and the energy will all have different orders).
To understand how the energy behaves though we need to move to higher orders.
For example for $\beta = 1$ and $n=1$ we have $$\phi(r) = ar + br^2 + cr^3 + o(r^3).$$ Plugging this into the equation we get $$-b - 3cr + \alpha (a + br) = E a r + o(r)$$ and so $b = \alpha a$ and $E = -3c/a + \alpha^2$.
For $\beta = 2$, $\alpha = -1$ and $n=2$ we have $$\phi(r) = ar^2 + br^3 + cr^4 + o(r^4)$$ and consequently $$a + 3br + 12cr^2 - (a + br + cr^2) = Ear^2 + o(r^2)$$ which means that $b = 0$ and $E = 11c/a$.
If we further assumed that $c/a$ is negative this would yields an answer to your question. I have no idea why this should be true unfortunately. Can someone confirm or refute this?