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$:

${\int }_{-1}^1{\int }_{-\sqrt{1-x^2}}^{{\sqrt{1-x^2}}_{}}{\int }_{x^2\+y^2}^1\left(y^3\+z\+x^2\right) \mathit{\ⅆ}z \mathit{\ⅆ}y \mathit{\ⅆ}x$ = $\frac{5}{12} \mathrm{\π}$ 

To compute the flux through $R$, note that there are two boundaries, the paraboloid, and the unit disk in the plane $z\=1$. To compute the flux through the paraboloid, note that on the paraboloid

If this be integrated over the unit disk, the result is

${\int }_{-1}^1{\int }_{-\sqrt{1-x^2}}^{{\sqrt{1-x^2}}_{}}\left(2 x^2{y}^3\+\left(2 y^2-x^2\right) \left(x^2\+y^2\right)\right) \mathrm{dy} \mathrm{dx}\=\frac{\mathrm{\π}}{6}$ 

On the upper boundary (disk), the outward normal is $\mathbf{N}\=\mathbf{k}$, so $\mathbf{F}·\mathbf{N}\=x^{2}z$, which becomes $x^2$ in the plane $z\=1$. Implementing the flux integral in polar coordinates gives

${\int }_0^{2 \mathrm{\π}}{\int }_0^{1}r\cdot {\left(r \cos \left(\mathrm{\θ}\right)\right)}^2 \mathrm{dr} d\mathrm{\θ}\=\frac{\mathrm{\π}}{4}$ 

The total flux is then $\left(\frac{1}{6}\+\frac{1}{4}\right) \mathrm{\π}\=\frac{5}{12} \mathrm{\π}$, the same value obtained for the volume integral of the divergence, as predicted by the Divergence theorem.

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 paraboloid and a unit disk in the plane $z\=1$. For the flux through the paraboloid, use a task template.

For the "open" surface $z\=x^2\+y^2$, Maple uses the normal ${\mathbf{R}}_x\times {\mathbf{R}}_y$, where $\mathbf{R}\=x \mathbf{i}\+y \mathbf{j}\+\left(x^2\+y^2\right) \mathbf{k}$ is a position-vector description of the paraboloid. This gives a normal that is upward, and hence inward for the closed region $R$. Because $R$ is a closed region, the normal should be outward, and hence downward; the flux through the paraboloid is actually $\mathrm{\π}\/6$.$$

On the upper boundary (disk), the outward normal is $\mathbf{N}\=\mathbf{k}$, so $\mathbf{F}·\mathbf{N}\=x^{2}z$, which becomes $x^2$ in the plane $z\=1$. The flux through this disk is given by the following calculation in which the integration is implemented in polar coordinates.

The total flux is then $\left(\frac{1}{6}\+\frac{1}{4}\right) \mathrm{\π}\=\frac{5}{12} \mathrm{\π}$, the same value obtained for the volume integral of the divergence, as predicted by the Divergence theorem.

For the "open" surface $z\=x^2\+y^2$, Maple uses the normal ${\mathbf{R}}_x\times {\mathbf{R}}_y$, where $\mathbf{R}\=x \mathbf{i}\+y \mathbf{j}\+\left(x^2\+y^2\right) \mathbf{k}$ is a position-vector description of the paraboloid. This gives a normal that is upward, and hence inward for the closed region $R$. Because $R$ is a closed region, the normal should be outward, and hence downward; the flux through the paraboloid is actually $\mathrm{\π}\/6$.$$

On the upper boundary (disk), the outward normal is $\mathbf{N}\=\mathbf{k}$, so $\mathbf{F}·\mathbf{N}\=x^{2}z$, which becomes $x^2$ in the plane $z\=1$. The integration of $x^2$ over the unit disk is simple enough to implement directly; the Flux command is a viable alternative.

The total flux is then $\left(\frac{1}{6}\+\frac{1}{4}\right) \mathrm{\π}\=\frac{5}{12} \mathrm{\π}$, the same value obtained for the volume integral of the divergence, as predicted by the Divergence theorem.