Q:

Why are the radial wavefunction and radial distribution function different?

I'm having trouble understanding the following two graphs, which are the radial wavefunction of the hydrogen 1s orbital and the corresponding radial distribution function:

Specifically, why is the radial distribution function zero at $r = 0$, but not the radial wavefunction?
Also, why does the radial distribution function have a maximum at the Bohr radius $a_0$?
(Mathematica input for the graphs: Plot[{2*Exp[-x], 4*Exp[-2 x]*x^2}, {x, 0, 6}, AxesLabel -> {r/Subscript[a, 0], Subscript[f, "1s"][r]}, PlotRange -> {{0, 6}, {0, 2}}, PlotLegends -> {"Radial wavefunction", "Radial distribution function"}, BaseStyle -> {FontSize -> 14}, ImageSize -> Medium])

A:

The radial distribution has a different form due to integration over the angles:
If we take the absolute square of the wave function $\Psi$ and integrate over the whole volume, we get the Norm of the wave function
\begin{equation}
N = \int \limits _{\varphi =0}^{2\pi }\ \int \limits _{\theta =0}^{\pi }\ \int \limits _{r=0}^{\infty } |\Psi(r,\theta ,\varphi )|^2r^{2}\sin \theta \,\mathrm {d} r\,\mathrm {d} \theta \,\mathrm {d} \varphi
\end{equation}
where $r^{2}\sin \theta \,\mathrm {d} r\,\mathrm {d} \theta \,\mathrm {d} \varphi$ comes from transforming the infinitesimal volume element from cartesian coordinates ($\mathrm{d} V = \mathrm {d} x\,\mathrm {d} y\,\mathrm {d} z$) to spherical coordinates.
For the radial probability distribution $P(r)$ we do the same, but omit the integral over $r$
\begin{equation}
P(r) = \int \limits _{\varphi =0}^{2\pi }\ \int \limits _{\theta =0}^{\pi }\ |\Psi(r,\theta ,\varphi )|^2r^{2}\sin \theta \,\mathrm {d} \theta \,\mathrm {d} \varphi
\end{equation}
which we can rearrange to 
\begin{equation}
P(r) = |R(r)|^2r^{2}\int \limits _{\theta =0}^{\pi }\int \limits _{\varphi =0}^{2\pi }|Y(\theta ,\varphi)|^2\sin \theta \,\mathrm {d} \theta \,\mathrm {d} \varphi
\end{equation}
Where $R(r)$ is the radial part of the wave function $\Psi(r,\theta ,\varphi)=R(r)Y(\theta ,\varphi)$.
For the $1s$ orbital we $Y(\theta ,\varphi)=1$ and straight forward integration over $\theta$ and $\varphi$ yields 
\begin{equation}
P(r) = |R(r)|^2r^{2} 4\pi
\end{equation}
So no matter what orbital $\Psi$ represents, $P(r)$ always has a factor $r^2$, which means we always have $P(r=0)=0$.
Your given example is the $1s$ orbital $\Psi_{1s}\propto \exp(-r)$. In contrast to all other orbitals, it does not have a polynomial in $r$. Therefore it has no node at $r=0$.

A:

The radial wave function $R(r)$ is simply the value of the wave function at some radius $r$, and its square is the probability of the finding an electron in some infinitesimal volume element around a point at distance $r$ from the nucleus.
But, the infinitesimal volume of space at radius $r$ is $4\pi r^{2} dr$ (it's a spherical shell with thickness $dr$ at radius $r$). That means that probability of finding an electron at radius $r$ is proportional to $R^{2}(r)4\pi r^{2}$. But the behavior of this function is such that the probability of finding the electron at radius 0 is also 0.

