The divergence of F:

$\nabla ·\mathbf{F}\={\partial }_x\left(x y^3\right)\+{\partial }_y\left(y z\right)\+{\partial }_z\left(x^{2}z\right)\=y^3\+z\+x^2$ 

Implement the integral of $\nabla ·\mathbf{F}$ over the interior of $R$ in spherical coordinates:

${\int }_0^{2 \mathrm{\pi }}{\int }_0^{\mathrm{\π}\/2}{\int }_0^1\left({\left(\mathrm{\ρ} \sin \left(\mathrm{\φ}\right)\sin \left(\mathrm{\θ}\right)\right)}^3\+\left(\mathrm{\ρ} \cos \left(\mathrm{\φ}\right)\right)\+{\left(\mathrm{\ρ} \sin \left(\mathrm{\φ}\right)\cos \left(\mathrm{\θ}\right)\right)}^2\right) {\mathrm{\ρ}}^2\sin \left(\mathrm{\φ}\right) d\mathrm{\ρ} d\mathrm{\φ} d\mathrm{\θ}$ = $\frac{23}{60} \mathrm{\π}$ 

To compute the flux through $R$, note that there are two boundaries, the upper hemisphere, and the disk that is its projection onto the plane $z\=0$. To compute the flux through the upper hemisphere, note that on that surface

$\mathbf{F}·\mathbf{N} \mathrm{d\σ}\=\left[\begin{array}{c}x y^3\\ y z\\ x^{2}z\end{array}\right]·\left[\begin{array}{c}x\\ y\\ z\end{array}\right] \frac{\mathrm{dA}}{z}$ = $\left(\frac{x^2{y}^3}{\sqrt{1-x^2-y^2}}\+y^2\+x^2\sqrt{1-x^2-y^2}\right) \mathrm{dA}$ 

If this be integrated over the unit disk in polar coordinates, the result is

${\int }_0^1{\int }_0^{2 \mathrm{\π}}r \left(\frac{r^5{\cos }^2\left(\mathrm{\theta }\right){\sin }^3\left(\mathrm{\theta }\right)}{\sqrt{1-r^2}}\+r^2{\sin }^2\left(\mathrm{\theta }\right)\+r^2{\cos }^2\left(\mathrm{\theta }\right)\sqrt{1-r^2}\right) d\mathrm{\θ} \mathrm{dr}$ = $\frac{23}{60} \mathrm{\π}$

On the lower boundary (disk), the outward normal is $\mathbf{N}\=-\mathbf{k}$, so $\mathbf{F}·\mathbf{N}\=-x^{2}z$, which becomes zero in the plane $z\=0$. Hence, the flux through the disk at the bottom of $R$ vanishes, and the flux through the upper surface matches the volume integral of the divergence.

The Student VectorCalculus package is needed for calculating the divergence, but it then conflicts with any multidimensional integral set from the Calculus palette. Hence, the Student MultivariateCalculus package is installed to gain Context Panel access to the MultiInt command.

There are two parts to the boundary of $R$, the surface of the upper hemisphere, and, in the plane $z\=0$, the unit disk that is the projection of the upper hemisphere. For the flux through the upper surface, use a task template.

On the lower boundary (disk), the outward normal is $\mathbf{N}\=-\mathbf{k}$, so $\mathbf{F}·\mathbf{N}\=-x^{2}z$, which becomes zero in the plane $z\=0$. Hence, the flux through the disk at the bottom of $R$ vanishes, and the flux through the upper surface matches the volume integral of the divergence.

On the lower boundary (disk), the outward normal is $\mathbf{N}\=-\mathbf{k}$, so $\mathbf{F}·\mathbf{N}\=-x^{2}z$, which becomes zero in the plane $z\=0$. Hence, the flux through the disk at the bottom of $R$ vanishes, and the flux through the upper surface matches the volume integral of the divergence.

Maple takes a noticeable amount of time to execute the integrations in Cartesian coordinates. In the Interactive section, these integrations are performed in spherical coordinates, and the time taken is insignificant.