The solid harmonics were explained by Misha, so let me fill in the rest of the details.
Laplacian operator is given by $\Delta = \partial_x^2 + \partial_y^2 + \partial_z^2$. First suppose there we are only interested in the radial part. Using chain rule (and letting $D \equiv \partial_r$, as in your references), we can write
$$\partial_x = r_{,x} D, \quad \partial_x^2 = r_{,xx} D + (r_{,x} D)^2.$$
We need to compute
$$ r_{,x} = {\partial \sqrt{x^2 + y^2 + z^2} \over \partial x } = {x \over r} $$
and
$$ r_{,xx} = \partial_x {x \over r} = {1 \over r} - {x^2 \over r^3}.$$
The expressions for $\partial_y$ and $\partial_z$ are of course similar, so that
$$\Delta = ({3 \over r} - {r^2 \over r^3})D + D^2 = {2 \over r}D + D^2.$$
Now, the spherical part of the Laplacian operator is given by $-{L^2 \over r^2}$ where $\mathbf L$ is the angular momentum operator. If we use spherical harmonics of level $l$ (which are eigenstates of L^2 corresponding to eigenvalue $l(l+1)$) and make a substitution $\phi \to {\phi \over r}$ we get the Griffith's result.