Check the divergence theorem for the function
\[\mathbf{v}=r^2 \cos \left( \theta \right) \hat{\mathbf{r}} +r^2 \cos \left( \phi \right) \hat{\mathbf{\theta }} -r^2 \cos \left( \theta \right) \sin \left( \phi \right) \hat{\mathbf{\phi }} \]Using as your volume one octant of the sphere of radius $R $. Make sure you include the entire surface.
Answer $\pi R^4 /4$
Let's start by recalling the divergence theorem (also known as The Fundamental Theorem for Divergences):
\[\int _{\mathcal{V}} \left( \nabla \cdot \mathbf{v} \right) d\tau =\oint _{\mathcal{S}} \mathbf{v}\cdot d\mathbf{a} \]In other words, if we want to integrate the divergence of a field over a volume, we could equally compute the (closed) surface integral of the field, where the closed surface is the surface that encloses the aforementioned volume.
We will integrate over an octant of a sphere of radius $R $. Feel free to move around the following 3D render of the octant to get familiar with its geometry.
To begin with the computations, let's take the divergence of $\mathbf{v} $:
\begin{align*} \nabla \cdot \mathbf{v} &= \frac{1}{r^2 }\frac{\partial }{\partial r}(r^2 v_r)+\frac{1}{r \sin \left( \theta \right) }\frac{\partial }{\partial \theta }(\sin \left( \theta \right)v_{\theta } )+\frac{1}{r\sin \left( \theta \right) }\frac{\partial v_{\phi }}{\partial \phi } \\ &= \frac{1}{r^2 }\frac{\partial }{\partial r}(r^2 r^2 \cos \left( \theta \right) )+\frac{1}{r \sin \left( \theta \right) }\frac{\partial }{\partial \theta }(\sin \left( \theta \right)r^2 \cos \left( \phi \right) )+\frac{1}{r\sin \left( \theta \right) }\frac{\partial }{\partial \phi } \left( -r^2 \cos \left( \theta \right) \sin \left( \phi \right) \right) \\ &= \frac{1}{r^2 } \left( 4r^3 \cos \left( \theta \right) \right) +\frac{1}{r \sin \left( \theta \right) } \left( r^2 \cos \left( \phi \right) \cos \left( \theta \right) \right) +\frac{1}{r \sin \left( \theta \right) } \left( - r^2 \cos \left( \theta \right) \cos \left( \phi \right) \right) \\ &= 4r\cos \left( \theta \right) +r \cot (\theta )\cos \left( \phi \right) -r \cot (\theta ) \cos \left( \phi \right) \\ &= 4 r \cos (\theta ) \end{align*}The volume integral is then (recall that $d\tau = r^2 \sin \left( \theta \right) drd\theta d\phi $):
\begin{align*} \int _{\mathcal{V}} (\nabla \cdot \mathbf{v})d\tau &= \int_{\phi =0}^{\pi /2} \int_{\theta =0}^{\pi /2} \int_{r=0}^{R} \left( 4r \cos \left( \theta \right) \right) r^2 \sin \left( \theta \right) drd\theta d\phi \\ &= \frac{\pi R^4 }{4} \end{align*}If you are not sure how to do the $\theta $-integral, consider using $u $-substitution with $u=\sin \left( \theta \right) $, which will turn your integral into $\int udu $.
Now, we have 4 surfaces to take care of to compute the surface integral.
Let's start by the surface at a constant radius (blue surface in the 3D render above). For it, we have that $d\mathbf{a}=r^2 \sin \left( \theta \right) d\theta d\phi \hat{\mathbf{r}} $. Thus this first surface integral has the value:
\begin{align*} S_1 &= \int_{\phi =0}^{\pi /2} \int_{\theta =0}^{\pi /2} \mathbf{v}\cdot (R^2 \sin \left( \theta \right) d\theta d\phi \hat{\mathbf{r}}) \\ &= \int_{\phi =0}^{\pi /2} \int_{\theta =0}^{\pi /2} \left( R^2 \cos \left( \theta \right) \right) (R^2 \sin \left( \theta \right) ) d\theta d\phi\\ &= \int_{\phi =0}^{\pi /2} \int_{\theta =0}^{\pi /2} R^4 \cos \left( \theta \right) \sin \left( \theta \right) d\theta d\phi\\ &= R^4 \int_{\phi =0}^{\pi /2} \int_{\theta =0}^{\pi /2} \cos \left( \theta \right) \sin \left( \theta \right) d\theta d\phi \\ &= R^4 \frac{1}{2}\int_{\phi =0}^{\pi /2} d\phi \\ &= \frac{R^4 \pi }{4} \end{align*}For the surface on the $xy $-plane, we have $d\mathbf{a}=rdrd\phi \hat{\theta } $, so:
\begin{align*} S_2 &= \int_{\phi =0}^{\pi /2} \int_{r =0}^{R} \mathbf{v}\cdot (r dr d\phi \hat{\theta }) \\ &= \int_{\phi =0}^{\pi /2} \int_{r =0}^{R} r^3 \cos \left( \phi \right) drd\phi \\ &= \frac{R^4 }{4} \end{align*}We are left with two surfaces at constant $\phi $, one at $\phi =0 $ and the other at $\phi =\pi /2 $. In both cases, we have $d\mathbf{a}=rdrd\theta \hat{\phi } $. Let's do the one at $\phi =0 $ first: \begin{align*} S_3 &= \int_{\theta =0}^{\pi /2} \int_{r =0}^{R} \mathbf{v}\cdot (r dr d\theta \hat{\phi }) \\ &= \int_{\theta =0}^{\pi /2} \int_{r =0}^{R} -r^3 \cos \left( \theta \right) \sin \left( 0 \right) drd\theta \\ &= \int_{\theta =0}^{\pi /2} \int_{r =0}^{R} 0 drd\theta \\ &= 0 \end{align*}
And the other one:
\begin{align*} S_4 &= \int_{\theta =0}^{\pi /2} \int_{r =0}^{R} \mathbf{v}\cdot (r dr d\theta \hat{\phi }) \\ &= \int_{\theta =0}^{\pi /2} \int_{r =0}^{R} -r^3 \cos \left( \theta \right) \sin \left( \pi /2 \right) drd\theta \\ &= \int_{\theta =0}^{\pi /2} \int_{r =0}^{R} -r^3 \cos \left( \theta \right) drd\theta \\ &= -\frac{R^4 }{4} \end{align*}Now we can finally write:
\begin{align*} \oint _{\mathcal{S}} \mathbf{v}\cdot d\mathbf{a} &= S_1 +S_2 +S_3 +S_4 \\ &= \frac{R^4 \pi }{4}+\frac{R^4 }{4}-\frac{R^4 }{4}+0\\ &= \frac{R^4 \pi }{4} \end{align*}Which agrees with the value obtained by integrating the divergence over the octant's volume.