Moment of Non-Uniform Sphere

$x = r \, cos \theta \, sin \phi \\ y = r \, sin \theta \, sin \phi \\ z = r \, cos \phi \\ 0 \leq r \leq 1 \\ 0 \leq \theta \leq 2 \pi \\ 0 \leq \phi \leq \pi$

Volume element:

$dV = r^2 \, sin \phi \, dr \, d \theta \, d \phi$

Differential mass:

$dm = \rho \, dV = r^3 \, \theta \, sin \phi \, dr \, d \theta \, d \phi$

Distance from z-axis:

$R^2 = x^2 + y^2 = r^2 \, sin^2 \phi$

Contribution to moment:

$dI = dm \, R^2 = r^5 \, \theta \, sin^3 \phi \, dr \, d \theta \, d \phi$

Total moment:

$I = \int_0^1 r^5 \, dr \, \int_0^{2 \pi} \theta \, d \theta \, \int_0^{\pi} sin^3 \phi \, d \phi \\ = \frac{1}{6} \, 2 \pi^2 \, \frac{4}{3} \\ = \frac{4}{9} \pi^2$

$\text{matrix}=\text{FullSimplify}\left[\text{Block}\left[\{x,y,z\},x=r \cos (\theta ) \sin (\phi );y=r \sin (\theta ) \sin (\phi );z=r \cos (\phi );\left( \begin{array}{ccc} y^2+z^2 & -x y & -x z \\ x (-y) & x^2+z^2 & -y z \\ x (-z) & y (-z) & x^2+y^2 \\ \end{array} \right)\right]\right] \Rightarrow \\ \left( \begin{array}{ccc} r^2 \left(\cos ^2(\phi )+\sin ^2(\theta ) \sin ^2(\phi )\right) & -r^2 \cos (\theta ) \sin (\theta ) \sin ^2(\phi ) & -r^2 \cos (\theta ) \cos (\phi ) \sin (\phi ) \\ -r^2 \cos (\theta ) \sin (\theta ) \sin ^2(\phi ) & r^2 \left(\cos ^2(\phi )+\cos ^2(\theta ) \sin ^2(\phi )\right) & -r^2 \cos (\phi ) \sin (\theta ) \sin (\phi ) \\ -r^2 \cos (\theta ) \cos (\phi ) \sin (\phi ) & -r^2 \cos (\phi ) \sin (\theta ) \sin (\phi ) & r^2 \sin ^2(\phi ) \\ \end{array} \right)$

$\text{jacobian}=\text{CoordinateTransformData}[\text{Spherical}\to \text{Cartesian},\text{MappingJacobianDeterminant}][\{r,\phi ,\theta \}]$

$\text{inertiaMatrix}=\int_0^{\pi } \left(\int_0^{2 \pi } \left(\int_0^1 \text{jacobian}\ \text{matrix}\ (\theta\ r) \, dr\right) \, d\theta \right) \, d\phi \Rightarrow \\ \left( \begin{array}{ccc} \frac{4 \pi ^2}{9} & \frac{\pi }{9} & 0 \\ \frac{\pi }{9} & \frac{4 \pi ^2}{9} & 0 \\ 0 & 0 & \frac{4 \pi ^2}{9} \\ \end{array} \right)$

$\text{Eigensystem}[\text{inertiaMatrix}] \Rightarrow \\ \left( \begin{array}{ccc} \frac{1}{9} \pi (1+4 \pi ) & \frac{4 \pi ^2}{9} & \frac{1}{9} \pi (-1+4 \pi ) \\ \{1,1,0\} & \{0,0,1\} & \{-1,1,0\} \\ \end{array} \right)$

Since the $z$ axis is one of the principal axes of the object, it was easy to read out the moment of inertia.