Technical detail in the solution of the hydrogen atom I'm trying to do an exercise in which you solve the Schrödinger equation for the hydrogen atom. Through the exercise, I've already shown that the wavefunction is:
$$ \psi_{n\ell m}(r,\theta,\varphi) = R_{n\ell}(r)Y^m_\ell (\theta,\varphi)$$
and that $Y^m_\ell (\theta,\varphi)$ are the spherical harmonics. Then, when solving for the radial part, the exercise tells me we need to study the asymptotic behaviour of $R(r)$ for large and small $r$. I had no problem showing that $R(r) \overset{\underset{+\infty}{}}{\sim} e^{-kr}$ but then the exercise gets completely bananas for the asymptotic part for small $r$. 
It tells me to introduce $u(r) \doteqdot rR(r)$, and then I managed to show the Schrödinger equation becomes:
$$-\frac{\hbar^2}{2m}{\mathrm{d}^2 u(r) \over \mathrm{d}r^2} - \frac{e^2}{4\pi\epsilon_0r}u(r)+\frac{\hbar^2}{2m}\frac{\ell(\ell+1)}{r^2}u(r) = E_n u(r)$$
Ok. Then the exercise tells me to suppose $u(r) \overset{\underset{0}{}}{\sim} r^\lambda$ and to prove $\lambda > -\tfrac{1}{2}$. No problem, this follows from the fact that $R(r)$ must be a $L^2$ function. Then, as long as $\ell \neq 0$, we can solve the above equation for $r\to 0$, showing that $\lambda = \ell + 1$. The next part of the exercise says:

Apply Stokes's theorem on a sphere to show that
$$\nabla^2 \left ( \frac{1}{r} \right ) = -4\pi\delta({\mathbf r}) $$
then use this result to prove that $\lambda = \ell + 1$ even when
  $\ell = 0$.

And I get completely stumped. I can prove the above formula, but I have no idea what this Laplacian has to do with anything.
 A: First, it is not true that $R(r)$ has to be $L^2$. Because the integration measure is $dV = r^2\cdot dr\cdot d\Omega$, and we integrate $|R|^2\cdot |Y|^2$ with this measure, it is $R(r)r=u(r)$ and not $R(r)$ itself that must be $L^2$.
Now, this is not just a correction of an unrelated minor mistake in your comments; it actually answers your main question. Why? Because the $L^2$ integrability would actually allow $u(r)$ to go like $r^0$ i.e. constant near $r\to 0$: the constant would clearly be square-integrable. However, if $u(r)$ were a constant, $R(r)$ would scale as $1/r$.
The whole wave function would go like $1/r$ near $r\to 0$; it would have this factor. Note that assuming $l=0$ which is the only case for which this singular $1/r$ behavior is imaginable, the angular part $Y$ goes like a constant, too. It is a completely legitimate candidate wave function, $\psi(r)\sim 1/r$. 
However, the wave function doesn't solve the Schrödinger's equation exactly because the Laplacian term in the equation produces a new term proportional to the delta-function (that's why the identity, justified by Stokes' theorem, is mentioned) which isn't canceled, and the equation therefore fails to hold. This wave function of the sort $1/r$ should better be checked separately because its Laplacian is "almost zero". However, when we check it carefully, we confirm that the Ansatz $R(r)\sim r^l$ is completely general, whether $l$ is zero or positive, and the $1/r$ wave function that obeys the required normalization conditions and that could be an "exceptional extra solution" actually isn't a solution.
A: This note is just a restatement of Lubos' answer (a tyro's attempt to understand his answer) and comes courtesy of Dirac:
To investigate the appropriate boundary conditions at $r=0$, consider the (different) situation of a free particle with $l=0$ and $E=0$.  In this case Schrodinger's equation reduces to Laplace's equation $\nabla^2 \Psi=0$, and the reduced one-dimension equation becomes $d^2u(r)/dr^2=0$.  
Now, while $R=1/r$ (or $u=rR=1$) obviously satisfies the 1-D equation, it does not satisfy the 3-D Laplace equation because, as you showed, a delta-function appears at the origin, violating Laplace's equation there. 
The point is that not every solution of the simple 1-D equation satisfies the 3-D Laplace's equation, and likewise, since the Laplacian is a term in the Schrodinger equation, not every solution of the "1-D Schrodinger equation" will satisfy the 3-D Schrodinger's equation.  In particular, as $r \rightarrow 0$, $R$ must not tend to $\infty$  as fast as $1/r$, or, equivalently, $rR=u \rightarrow 0$, in order to avoid the delta-function.
Not quite done, since we haven't showed $u$ goes like $r$ near $r=0$, just that it goes to 0.  Assume a solution with the proper large-$r$ behavior ($u= \sum c_sr^s \exp{-kr}$), where consecutive values of $s$ differ by unity, although these values may not be integers.  From the above, the minimum $s$ must be positive.  Sub'ing and plugging, one finds that the minimum of $s$ must be either $l+1$ or $-l$. The lower value is ruled out, even when $l=0$, since it is not positive, so u goes like $r^{l+1}$ in all cases, including $l=0$.
I'll reiterate that I've just restated Lubos' answer (which I now understand!).
