Solution (source code)

= Solution

For $i\ne j$, the <inertia tensor> component is $I_{ij}=-\rho\int_Vx_ix_j\,dV$. Reflect either coordinate appearing once: the cylindrical domain and density are unchanged, while the integrand changes sign. Thus every off-diagonal component vanishes. Rotational symmetry also implies $I_{11}=I_{22}$.

Use <cylindrical coordinates> $x_1=r\cos\theta$, $x_2=r\sin\theta$, $x_3=z$, with volume element $r\,dr\,d\theta\,dz$. The axial <moment of inertia> is
$$
I_{33}=\rho\int_{-2}^{2}dz\int_0^{2\pi}d\theta\int_0^1r^3\,dr=2\pi\rho.
$$
The two contributions to a transverse component are
$$
\int_Vx_2^2\,dV=4\left(\int_0^1r^3dr\right)\left(\int_0^{2\pi}\sin^2\theta\,d\theta\right)=\pi,
\qquad\int_Vz^2\,dV=\pi\int_{-2}^{2}z^2dz=\frac{16\pi}{3}.
$$
Therefore
$$
\boxed{I=\pi\rho\,\operatorname{diag}\left(\frac{19}{3},\frac{19}{3},2\right)=M\operatorname{diag}\left(\frac{19}{12},\frac{19}{12},\frac12\right),\qquad M=4\pi\rho.}
$$
The axial value $MR^2/2$ and transverse value $M(3R^2+h^2)/12$ for $R=1$, $h=4$ give an independent check on both diagonal entries.