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 }_{-7\/4}^{7\/4}{\int }_{y_{-}}^{y_{\+}}{\int }_{x^2\+4 y^2}^{3-y}\left(y^3\+z\+x^2\right) \mathrm{dz} \mathrm{dy} \mathrm{dx}$ = $\frac{18917479}{3145728}\mathrm{\π}$

To compute the flux through $R$, note that there are two boundaries, the upper one being the plane $z\=3-y$, and the lower one being the elliptic paraboloid $z\=x^2\+4 y^2$. To compute the flux through the lower surface, note that on that surface

where N points downward on the paraboloid, and thus outward with respect to $R$.

If this be integrated over the ellipse $x^2\+4 y^2\+y\=3$, the result is

${\int }_{-7\/4}^{7\/4}{\int }_{y_{-}}^{y_{\+}}\left(2 x^2{y}^3\+\left(8 y^2-x^2\right) \left(x^2\+4 y^2\right)\right) \mathrm{dy} \mathrm{dx}$ = $\frac{10196263}{3145728}\mathrm{\π}$ 

On the upper boundary (the plane), the upward (and hence outward) normal is $\mathbf{N}\=\left(\mathbf{j}\+\mathbf{k}\right)\/\sqrt{2}$, so

$\mathbf{F}·\mathbf{N} \mathrm{d\σ}\=\left(\left(x^2\+y\right)z\/\sqrt{2}\right) \left(\sqrt{2}\mathrm{dA}\right)\=\left(x^2\+y\right)\left(3-y\right) \mathrm{dA}$ 

on that surface. If this be integrated over the ellipse $x^2\+4 y^2\+y\=3$, the result is

${\int }_{-7\/4}^{7\/4}{\int }_{y_{-}}^{y_{\+}}\left(x^2\+y\right)\left(3-y\right) \mathrm{dy} \mathrm{dx}$ = $\frac{45423}{16384}\mathrm{\pi }$

The total flux through the boundaries of the region $R$ is then the sum

$\frac{45423}{16384}\mathrm{\π}\+\frac{10196263}{3145728}\mathrm{\π}$ = $\frac{18917479}{3145728}\mathrm{\π}$ 

which matches the volume integral of the divergence of F inside $R$.

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 plane $z\=3-y$, and the elliptic paraboloid $z\=x^3\+4 y^2$. For the flux through the upper surface, use a task template.

On the lower bounding surface, namely, the elliptic paraboloid, the flux can be computed by the same task template. However, because this is an "open" surface, Maple will default to the normal ${\mathbf{R}}_x\times {\mathbf{R}}_y$, where R is a position-vector representation of the paraboloid. This normal will be upward on the surface, but inward with respect to the region $R$. Hence, the sign of the computed flux will have to be changed to get the flux in the direction of the outward normal.

Making the appropriate sign change in the flux through the paraboloid, the total flux through the boundaries of the region $R$ is then the sum

$\frac{45423}{16384}\mathrm{\π}\+\frac{10196263}{3145728}\mathrm{\π}$ = $\frac{18917479}{3145728}\mathrm{\π}$ 

which matches the volume integral of the divergence of F inside $R$.

Because the lower boundary is an "open" surface, Maple defaults to the normal ${\mathbf{R}}_x\times {\mathbf{R}}_y$, where R is a position-vector representation of the paraboloid. This normal will be upward on the surface, but inward with respect to the region $R$. Hence, the sign of the computed flux will have to be changed to get the flux in the direction of the outward normal.

The flux through the upper boundary can also be obtained with the Flux command.

Making the appropriate sign change in the flux through the paraboloid, the total flux through the boundaries of the region $R$ is then the sum

$\frac{45423}{16384}\mathrm{\π}\+\frac{10196263}{3145728}\mathrm{\π}$ = $\frac{18917479}{3145728}\mathrm{\π}$ 

which matches the volume integral of the divergence of F inside $R$.