Solution (source code)

= Solution

Take $0<r_1<r_2$, so the shell excludes the possible singularity at the origin. The <divergence and curl of a radial vector field> formula in <spherical coordinates> gives
$$
\nabla\cdot\mathbf F
=\frac1{r^2}\frac{d}{dr}\left(r^2r^{-\beta}\right)
=(2-\beta)r^{-\beta-1}.
$$
With <volume element> $dV=r^2\sin\theta\,dr\,d\theta\,d\phi$, direct integration gives
$$
\begin{aligned}
\iiint_V\nabla\cdot\mathbf F\,dV
&=\int_{r_1}^{r_2}\int_0^\pi\int_0^{2\pi}
(2-\beta)r^{1-\beta}\sin\theta\,d\phi\,d\theta\,dr\\
&=\boxed{4\pi\left(r_2^{2-\beta}-r_1^{2-\beta}\right).}
\end{aligned}
$$
For the independent boundary calculation, the outer <sphere> has outward <unit normal vector> $+\hat{\mathbf r}$ and the inner <sphere> has outward <unit normal vector> $-\hat{\mathbf r}$. Their <fluxes> are
$$
\Phi_{\rm outer}=4\pi r_2^2r_2^{-\beta},\qquad
\Phi_{\rm inner}=-4\pi r_1^2r_1^{-\beta}.
$$
Their sum equals the direct <volume integral>, verifying the <divergence theorem>. The restriction $\beta\ne2$ permits the displayed power-antiderivative step. At $\beta=2$, treated separately, both the shell <volume integral> and net boundary <flux> are zero; the excluded origin must still not be filled in.