Differentials in Spherical Shell - Maxwell Distribution In explaining the Maxwell distribution of molecular speeds, my pchem textbook uses the following figure: 
We are basically trying to find the probability of having a particle with a speed $u$ between $u$ and $u + du$. The axes are the speed components in each spatial dimension. In the book's explanation, they state that the volume of this shell, $$ 4 \pi u^2 du = du_x du_y du_z . $$ Could someone help explain why this is true? The only substitution I can think to make is: $$  4 \pi u^2 du = 4 \pi (u_x^2 + u_y^2 + u_z^2)\sqrt{du_x^2+du_y^2+du_z^2}, $$ which doesn't really get me anywhere.
 A: The intuition is that $4\pi u^2$ is the area of a sphere of radius $u$, and now to find the volume of the thin shell between radius $u$ and radius $u+du$, you multiply the area of the surface of the shell by the thickness of the shell and find that its volume is $4\pi u^2du$.
We can make this intuition more precise in a number of ways.  Here are two:
Method 1. (Algebraic)
Recall that the volume of a sphere of radius $u$
\begin{align}
  V(u) = \frac{4}{3}\pi u^3
\end{align}
Now, this means that for any $\Delta u>0$, the volume of a sphere of radius $u+\Delta u$ is
\begin{align}
  V(u+\Delta u) 
&= \frac{4}{3}\pi (u+\Delta u)^3 \\
&= \frac{4}{3}\pi(u^3 + 3u^2\Delta u+3u\Delta u^2 + \Delta u^3)
\end{align}
The difference between these volumes gives us the volume of the shell between radii $u$ and $u+\Delta u$.
\begin{align}
  V(u+\Delta u) - V(u) = 4\pi u^2\Delta u + \left( 4u\Delta u^2 + \frac{4}{3}\pi\Delta u^3\right)
\end{align}
but notice that for $\Delta u$ small and $u$ fixed, the first term on the right dominates the other two which are higher powers of $\Delta u$, so the expression $4\pi u^2\Delta u$ is the first order approximation to the shell volume.  This first order approximation is precisely the thing that counts in an integral.
Method 2. (Integration)
The volume of a given region $R$ in three dimensions can be written in cartesian coordinates as a multiple integral over that region as follows:
\begin{align}
  \mathrm{vol}(R) = \int_R du_x\,du_y\,du_z
\end{align}
Now, if you perform the change of variables to spherical coordinates, you find
\begin{align}
  u_x &= u\sin\theta\cos\phi \\
  u_y &= u\sin\theta\sin\phi \\
  u_z &= u\cos\theta
\end{align}
which, as you can check, gives
\begin{align}
  du_x\,du_y\,du_z = u^2\sin\theta\,du\,d\theta\,d\phi \tag{$\star$}
\end{align}
so the volume of the given region can be written as
\begin{align}
  \mathrm{vol}(R) = \int_R u^2\sin\theta\,du\,d\theta\,d\phi
\end{align}
Note that one needs to be careful to write the appropriate limits of integration in each case that correctly specify the region $R$ in the given system of coordinates.
Now, in spherical coordinates, a shell with inner and outer radii $u$ and $u+\Delta u$ is the region of points $(u, \theta, \phi)$ satisfying the addition constraints
\begin{align}
  0\leq\theta\leq\pi, \qquad 0\leq\phi<\pi/2
\end{align}
so the volume of this region is
\begin{align}
  \mathrm{vol}(\mathrm{shell}) 
&= \int_0^{2\pi}\int_0^\pi\int_u^{u+\Delta u} u^2\sin\theta\,du\,d\theta\,d\phi\\
&= \int_0^{2\pi} d\phi\int_0^\pi \sin\theta\,d\theta\int_u^{u+\Delta u} u^2du
\end{align}
Now recall that
\begin{align}
  \int_0^{2\pi} d\phi = 2\pi, \qquad \int_0^\pi \sin\theta\,d\theta = 2 
\end{align}
so that
\begin{align}
  \mathrm{vol}(\mathrm{shell}) = 4\pi \int_u^{u+\Delta u} u^2du = \frac{4}{3}\pi (u+\Delta u)^3 - \frac{4}{3}\pi u^3
\end{align}
But wait!  This is just what we got in the "algebraic" method above.  Of course, all we've really done here is derive the formula for the volume of the sphere using integration.  The real magic was in the change of variables which gave eq. $(\star)$
A: 
"What is the probability that a particle's speed lies between $u$ and $u+du$?"

Short answer
The probability that a particle's velocity has an $x$-component between $u_x$ and $u_x+du_x$, and similarly for $y$ and $z$, is proportional to the volume taken up by that small cubic chunk of velocity-space, $du_xdu_ydu_z.$ Given these three components of the velocity you can calculate the speed and direction (within a small range) that the particle is traveling.
Integrating over all possible directions,
$$\int_{\theta,\phi} du_xdu_ydu_z = \int_{\theta,\phi} u^2\sin(\theta)d\theta d\phi du=4\pi u^2du.$$
Details
The probability that the $x$-velocity lies somewhere between $a$ and $b$ is given by
$$P_{u_x}(a,b)=\int_{a}^{b}f(u_x)du_x$$
where $f(u_x)$ is called the probability distribution of $u_x$. By the continuity of the integral, we see that
$$P_{u_x}(u_x,u_x+dx)=\int_{u_x}^{u_x+dx}f(u_x)du_x=f(u_x)du_x=dP_{u_x}$$
Since the magnitude of the $x$-velocity cannot affect the magnitude of the $y$-velocity and so on, the probabilities multiply, and the probability that the particle has all three components in some range is
$$dP_{u_x}dP_{u_y}dP_{u_z}=f(u_x)f(u_y)f(u_z)du_x du_y du_z \tag{1}$$
If we want to count up all the velocity states with the same velocity exclusive of direction, then we write
\begin{align}
P_u(u,u+du) &= \int_{\theta,\phi} f(u_x)f(u_y)f(u_z) du_xdu_ydu_z \\
&= \int_{\theta,\phi} f(u) u^2\sin(\theta)d\theta d\phi du \\
&=f(u)4\pi u^2du.\tag{2}
\end{align}
Comparing equations (1) and (2) should hopefully remove the confusion.
Note: $f(u)$ is distinct from the $f(u_i)$, but this should be clear from the context.
