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.
I think this problem is similar to the problem of finding modes of rectangular dielectric waveguide. In this case, you can use the effective-index method for finding the approximated solution (For your problem, we can call it effective-potential method).
For more information about effective-index method see the following articles:
The basis of this method is that the mode of a waveguide can be separated into products of two functions, one in $x$ direction which is dependent only on $x$ and one in $y$ direction which is dependent only on $y$. These can be solved independently and combined to produce the mode structure. In this way, the 2D waveguide structure can be separated into two single structures, one being a step index planar waveguide in $x$ direction and other in $y$ direction. In fact, this is same as your suggestion for introducing $V_x$ and $V_y$, but in a special way that the solution is very closed to the exact solution
Best Answer
The most general widely studied examples of potentials are those that are derived from an exactly given ground state, and which have the property that they contain enough parameters to be closed under taking supersymmetric conjugates. These are called "shape invariant" potentials.
Given a real positive ground state, one can ask "Which potential has this ground state?" The answer is that if the ground state is $\exp(-W(x))$, and its energy is exactly zero, then the potential is:
$ V(x) = {1\over 2} |\nabla W|^2 - {1\over 2} \nabla^2 W$
The conjugate potential is defined with a plus sign between the two terms instead of a minus sign. It is more usual to define W as the derivative of what I am calling W, but this convention is terrible in higher than one dimension.
The two potentials taken together define a supersymmetric quantum mechanics, as originally defined by Witten. The supersymmetric quantum mechanics in imaginary time is a stochastic Brownian process with a drift which is an analytic function of the position. The conjugate potentials correspond to reversing the direction of the drift, and their properties are similar because they are related by a stochastic version of time-reversibility.
If W(x) goes to plus infinity at infinity (so that it actually defines a normalizable ground state), then the ground state is unique, and the conjugate potential has the exact same spectrum as the original potential, except it omits the lowest energy state. This, plus the form of the supercharge, gives exact solutions of many classes of quantum potentials.
Here are some simple W's which correspond to usual elementary quantum mechanics examples:
http://arxiv.org/abs/hep-th/9405029 has a bunch of more interesting examples. Any quantum mechanical potential which has closed form energy states is in this class.
A completely diffeent class of widely studied potentials are random potentials, as studied by Halperin and others, to understand Anderson localization.
Later Edit: The original paper by Anderson which started the random potential field is http://prola.aps.org/abstract/PR/v109/i5/p1492_1 "Absence of diffusion in certain random lattices", and it's one of the great classics. The setup is a square lattice with a random potential at each site, an independent random number between -V and V. The continuum limit in one dimension, where the potential is a random gaussian at each point is analyzed by Halperin B. I. Halperin, Green ' s Functions for a Particle in a One-Dimensional Random Potential, Phys. Rev. 139 , A104 (1965). The field is enormous--- look up "localization" on google scholar. It includes "weak localization" effects, which were popular in the mid 90s because they imply that resistance can drop sharply in the presence of a magnetic field, because the perturbative precurser to the localization process is hindered.