In multi-electron atoms the energy of the s, p and d orbitals for a given $n$ is determined primarily by shielding by inner electrons and penetration of the orbital towards the nucleus. As s orbitals have some probability of being at the nucleus their energy is lowest, next is p then d orbitals in increasing energy. Also there is spin-orbit coupling which splits sub-levels of p and d orbitals or any orbital which has angular momentum so this does not effect s orbitals. This effect is much smaller in d than p orbitals and generally smaller than the separation of s to p to d.
The statements in the book extract you give seem contradictory. The mean radial value of an H type orbital is calculated using $\displaystyle \int_0^\infty r^3R_{n,l}^*R_{n,l} dr $ where $R_{n,l}$ is the (normalised) radial part of the wavefunction for an atom with quantum numbers $n,\;l$. The $r^3$ arises because the volume element in the integration is $r^2\sin(\theta)d\phi d\theta$ in spherical coordinates. The other $r$ is how the average is defined, i.e. $\int r\psi^*\psi d\tau$ for some normalised wavefunction $\psi$.
The radial parts of the wavefunction for the 3s,3p and 3d are define with $q=Zr/a_0$ ($a_0$ is the Bohr radius) as
$$\begin{align}
R_{3,0}&=N(27-18q+2q^2)exp^{-q/3}\\
R_{3,1}&=N(6q-q^2)e^{-q/3}\\
R_{3,2}&=Nq^2e^{-q/3}
\end{align}$$
and $N$ is the normalisation which is different for each wavefunction. The integrals are not hard to evaluate but very tedious as there are many terms. However, there is a general formula for the average of H atom type orbitals and is
$$\langle r\rangle =\frac{n^2a_0}{Z}\left( 1+\frac{1}{2} \left( 1-\frac{l(l+1)}{n^2} \right) \right ) $$
The average values are
$$3s \;\langle r \rangle = 0.476 \text{ nm,}\quad 3p \;\langle r \rangle = 0.441\text{ nm,}\quad 3d\; \langle r \rangle = 0.370\text{ nm}$$
The s orbital value is the same as in your question but the others differ both in value and in order of size. Perhaps there are typos in the text you give or some other unstated factors producing these values.
The maximum values have to be calculated by making the derivative (of $q^2R_{n,l}^2$) zero and solving the equation. Alternatively the value can be taken with sufficient accuracy from a plot of the function which is far easier. The values are
$$3s \;r_{max} = 0.46 \text{ nm,}\quad 3p \;r_{max} = 0.43\text{ nm,}\quad 3d\; r_{max} = 0.32\text{ nm}$$
the plot below shows the radial functions $q^2R_{n,l}^2$ each normalised to the same maximum height. The mean is greater than the maximum in each case as expected as the main part of the radial probability is 'lopsided' towards larger distances.
(In evaluating the integrals $\displaystyle \int_0^\infty r^ne^{-\beta r}= \frac{n!}{\beta^{n+1}}$ is useful.)