The Divergence of the vector function comes out to be $4rcos(\theta$). The volume integral is thus: $$\int_{0}^R\int_{0}^{\pi/2}\int_{0}^{\pi/2} 4rcos(\theta)*(r^2sin(\theta) d\theta d\phi dr)=\pi R^4/4$$
Now, for finding the closed surface integral $\oint_{s}\vec{v}.d\vec{a}$, we write: $$\oint_{s}\vec{v}.d\vec{a}=\int_{1}\vec{v}.d\vec{a}+\int_{2}\vec{v}.d\vec{a}+\int_{3}\vec{v}.d\vec{a} +\int_{4}\vec{v}.d\vec{a} $$
Where $1$: the surface "bulging outward", $2,3,4$ are the surfaces (quarter-circles) in the $xy,yz,xz$ planes respectively.
Now, $\int_{1}\vec{v}.d\vec{a}=\int_{1}\vec{v}.(rd(\theta)rsin(\theta)d\phi\hat{r})=\int_{1}\vec{v}.(rd(\theta)rsin(\theta)d\phi\hat{r})=\int_{0}^{\pi/2}\int_{0}^{\pi/2}sin(\theta)cos(\theta)d\theta d\phi=\pi R^4/4$.
Which means, to "check" the divergence theorem , we need to show that $\int_{2}\vec{v}.d\vec{a}+\int_{3}\vec{v}.d\vec{a} +\int_{4}\vec{v}.d\vec{a} =0$.
However, I am finding this computation incredibly difficult, as the area vectors in these surfaces are pointing in the $\hat{k},\hat{j},\hat{i}$ directions, so we need to effectively convert between spherical and Cartesian coordinates.
Questions:
1. Is there a "symmetry"/alternate argument to show: $\int_{2}\vec{v}.d\vec{a}+\int_{3}\vec{v}.d\vec{a} +\int_{4}\vec{v}.d\vec{a} =0$?
2. Am I misinterpreting the $\oint$ symbol in this context? Is it somehow by definition only $\int_{1}$?
Based on your working, you are considering $\theta$ as the polar angle. I do not see it mentioned clearly in the question but if that is the notation used in the book and in the question, here is how the transformation would look -
$\begin{bmatrix} \hat{r} \\ \hat{\theta} \\ \hat{\phi} \end{bmatrix} = \begin{bmatrix} \sin\theta\cos\phi & \sin\theta\sin\phi & \cos\theta \\ \cos\theta\cos\phi & \cos\theta\sin\phi & -\sin\theta \\ -\sin\phi & \cos\phi & 0 \end{bmatrix} \begin{bmatrix} \hat{x} \\ \hat{y} \\ \hat{z} \end{bmatrix}$
So to show that the flux through $3$ planar parts of the sphere in the first octant sum to zero -
i) For part of the sphere in $YZ$ plane in first octant, the outward unit normal vector in cartesian coordinates will be $(-1, 0, 0)$ and as $\phi = \frac{\pi}{2}, 0 \leq \theta \leq \frac{\pi}{2} \,$ for the surface, it will transform to $(0, 0, 1)$ in spherical $(r, \theta, \phi)$. The vector field will be $(r^2 \cos\theta, 0, -r^2\cos\theta)$.
ii) For part of the sphere in $XZ$ plane in first octant, the outward unit normal vector in cartesian coordinates will be $(0, -1, 0)$ and as $\phi = 0, 0 \leq \theta \leq \frac{\pi}{2}$, it will transform to $(0, 0, -1)$ in spherical. The vector field will be $(r^2 \cos\theta, -r^2, 0)$. As the dot product is zero, the surface integral will be zero.
iii) For part of the sphere in $XY$ plane in first octant, the outward unit normal vector in cartesian coordinates will be $(0, 0, -1)$ and as $\theta = \frac{\pi}{2}, 0 \leq \phi \leq \frac{\pi}{2}$, it will transform to $(0, 1, 0)$ in spherical. The vector field will be $(0, r^2 \cos \phi, 0)$.
For $(i)$ and $(iii)$, you can either calculate and show or simply argue that the integral will be equal but with opposite sign and hence sum to zero.