First of all, this is not a homework question. The point of this exercise is to understand how Hydrogen atom, Positronium and ultimately strong bound states behave when Darwin-type potential is added to Coulomb potential. In all textbooks this is treated as perturbation theory problem (in the context of fine structure), but that fails even when perturbation theory converges, let alone when it doesn't. 

For example $U^{91+}$ atom has effective $\alpha=92/137\approx 2/3$ (so perturbation expansion should in principle converge), and the first order correction for Darwin potential is of order $\alpha^4$ which is about half the binding energy, about an order of magnitude larger than the actual value. One can calculate higher order contributions, but it becomes impossible to do it beyond 4th order, and even that doesn't really produce results (because many more terms are needed to get the desired accuracy). 

In this particular example one can use Dirac equation to get the energy, but that includes both relativistic corrections, spin-orbit corrections, Darwin corrections and it's impossible to separate the contributions. When $\alpha>1$, Dirac equation fails (formula for energies returns complex numbers) and that approach fails completely (also, perturbation expansion seems to diverge). 

So the idea is to first solve the Schroedinger equation for non-relativistic Hydrogen atom with delta-function perturbation analytically and get an understanding how that delta-function affects energies and wave functions when the interaction is strong, and hence separate effects of Darwin term from relativistic effects. This would in principle also answer (or if not answer, then shed some light on) the question why does the Dirac equation fail when $\alpha>1$? For $l=0$ solutions spin-orbit interaction doesn't contribute, so is the Darwin term or the relativistic corrections that that cause the problem? 

Therefore, I'm trying to solve the Schroedinger equation for Hydrogen atom with delta-function perturbation (essentially Darwin potential) 
$$
\frac{p^2}{2m}\psi -\frac{\alpha}{r}\psi +b \delta^3(\vec{r})\psi=E\psi
$$
where parameter $b$ is not necessarily small. If I express the 3D delta function in terms of 1D delta function
$$
\delta^3(\vec{r})=\frac{\delta(r)}{2\pi r^2}
$$
I get
$$
\frac{\partial^2 \psi}{\partial r^2} + \frac{2}{r}\frac{\partial \psi}{\partial r} +\frac{2m\alpha}{r}-\frac{m b/\pi\delta(r) +l(l+1)}{r^2}\psi +2mE\psi=0
$$
which is exactly the same as unperturbed Hydrogen atom for $r>0$, so the only difference is at $r=0$. From perturbation theory I know that the effect of $\delta$-function cannot be neglected, and is finite only for $l=0$ states (since unperturbed wave function vanishes for $l>0$ states at origin), but I'm uncertain how to solve this equation? 

For $s$ orbital solutions, effects of the $\delta$-function would presumably show as some sort of boundary condition at $r=0$? Since the added potential is infinitely positive at $r=0$, presumably wave function should vanish at $r=0$? Also, presumably there's some sort of condition on it's derivative? 

Sounds like I'm missing something here, but I don't know what.