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. It sounds like it should be reasonably easy, but all textbooks treat the $\delta$-term as a small perturbation that contributes to fine structure for $s$ states. 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.