When trying to solve the Schrödinger equation for hydrogen, one usually splits up the wave function into two parts:
$$\psi(r,\phi,\theta)= R(r)Y_{\ell,m}(\phi,\theta).$$
I understand that the radial part usually has a singularity for the $1s$ state at $r=0$ and this is why you remove it by writing:
$$R(r) = \frac{u(r)}{r}$$
But what is the physical meaning of
$$R(r=0) = \infty~?$$
Wouldn't this mean that the electron cloud is only at the centre of the atomic nucleus?
The infinitesimal probability for the electron to be in the volume $dV$ around a point $(r,\theta,\phi)\leftrightarrow (x,y,z)$ is given by $$ dP = dV\cdot |\psi(x,y,z)|^2 = dV\cdot |R(r)|^2\cdot |Y_{lm}(\theta,\phi)|^2 =\dots$$ as you can see if you substitute your Ansatz for the wave function. However, the infinitesimal volume $dV=dx\cdot dy\cdot dz$ may be rewritten in terms of differentials of the spherical coordinates as $$ dV = dr\cdot r^2 \cdot d\Omega = dr\cdot r^2 \cdot \sin\theta\cdot d\theta\cdot d\phi $$ where the small solid angle $d\Omega$ was rewritten in terms of the spherical coordinates. You see that for dimensional reasons (or because the surface of a sphere scales like $r^2$), there is an extra factor of $r^2$ in $dV$ and therefore also in $dP$ which suppresses the probability. There is simply not enough volume for small values of $r$.
So $|R(r)|^2$ may still go like $1/r^2$ for small $r$ and in that case, $dV$ will be proportional to $dr$ times a function that is finite for $r\to 0$. Such $dP$ may be integrated and there's no divergence at all near $r=0$.
That's why one should allow the wave function to go like $1/r$ near $r=0$ which is the true counterpart of one-dimensional wave function's being finite near a point. However, Nature doesn't use this particular loophole because the wave function $\psi$ for small $r$ actually scales like $r^l$ where $l$ is the orbital quantum number and the wave function actually never diverges even though it could.
Update 2016: I should have and could have written it four years ago but I didn't. While the normalizability allows $1/r$ around $r=0$, such singular functions ultimately can't be in stationary or nearly stationary states, for the following reason which differs from various reasons above and those in the comments.
For example, someone mentioned that $1/r$ could lead to a continuous spectrum or some surprising degeneracies. But if the correct wave functions predicted a continuous or degenerate spectrum in a box, then it would be how Nature works. The actual reason why $1/r$ is not finally allowed as a stationary wave function near $r=0$ is that the Laplacian of this wave function (and Schrödinger's equation contains such a Laplacian) is proportional to a delta-function at the origin (or contains such a term) and no other term in Schrödinger's equation can cancel this delta-function, so the Schrödinger equation must be violated.
The physical observable is not the wavefunction, but its integral over a finite area. In spherical coordinates, this is:
$$P({\vec x})=\int~\mathrm dr\, \mathrm d\theta\, \mathrm d\phi ~r^{2}\sin\theta~ \psi^{*}\psi$$
This integrand is manifestly finite at $r=0$, even if $R(r)$ has a $\frac{1}{r}$ divergence.
For a hydrogen-like atom in 3 spatial dimensions, the rewriting of the radial part
$$R(r)~=~\frac{u(r)}{r}$$
is not performed to keep the $u(r)$ part regular, as OP suggests, but usually because the 3D radial equation in terms of the $u$ function has the same form as a 1D Schrödinger equation.
Imagine that the radial wave function goes as a power
$$R(r) ~\sim ~ r^{p} \qquad {\rm for} \qquad r~\to~ 0, \qquad p~\in~\mathbb{R}.$$
On general grounds, one can impose the following list of consistency conditions, listed with the weakest condition first and the strongest condition last.
Normalizability of the wave function $$\infty~>~\langle\psi|\psi\rangle~=~\int \mathrm{d}^3r~|\psi(\vec{r})|^2 ~\propto~ \int_0^{\infty} r^{2}\mathrm{d}r~|R(r)|^2 .$$ Integrability at $r=0$ yields that the power $p>-\frac{3}{2}$. In other words, this normalizability condition does not by itself imply that $R(r)$ or $u(r)$ should be regular at $r=0$, which is also the conclusion of many of the other answers.
The expectation value of the potential energy $V$ should be bounded from below, $$-\infty~<~\langle\psi| V|\psi\rangle~=~\int \mathrm{d}^3r~V(r)|\psi(\vec{r})|^2~\propto~-\int_0^{\infty} r \, \mathrm{d}r~|R(r)|^2. $$ Integrability at $r=0$ yields that the power $p>-1$. In other words, $u(r)$ should be regular for $r\to 0$.
The kinetic energy operator (or equivalently, the Laplacian $\Delta$) should behave self-adjointly for two wave functions $\psi_1(\vec{r})$ and $\psi_2(\vec{r})$, $$\langle\psi_1| \Delta\psi_2\rangle~=~-\langle\vec{\nabla}\psi_1| \cdot\vec{\nabla}\psi_2\rangle,$$ without picking up pathological contributions at $r=0$. A detailed analysis shows that the powers of the radial parts of $\psi_1(\vec{r})$ and $\psi_2(\vec{r})$ should satisfy $p>-\frac{1}{2}$.
In comparison, the actual bound state solutions have non-negative $p=\ell\in \mathbb{N}_0$, and therefore satisfy these three conditions.
In addition to the simply geometric constraints that Jerry and Lubos talk about, the derivation used to illustrate the problem almost always assumes that the proton is a point particle which is a pretty good approximation but not strictly true. Working the problem again with a realistic proton charge density function (roughly constant inside a radius of about 1 fm) would be another way to remove the singularity.
Mind you, you this argument does not hold true for the positronium so you still need the geometric constraint.
For Hydrogen, $R(r)$ does not diverge, as $U(r)$ vanishes as fast as (or faster than) $r$ as $r\rightarrow 0$. In fact, it's only for the $s$ orbitals that the wavefunction is non zero at $r=0$. But as pointed out before, a non-zero radial wavefunction does not mean a non-zero probability of finding the electron at the center.
