Use the Divergence Theorem to evaluate the surface integral

Compute the divergence of [tex]\vec F[/tex].
[tex]\mathrm{div}(\vec F) = \dfrac{\partial(2x^3+y^3)}{\partial x} + \dfrac{\partial (y^3+z^3)}{\partial y} + \dfrac{\partial(3y^2z)}{\partial z} = 6x^2 + 3y^2 + 3y^2 = 6(x^2+y^2)[/tex]
By the divergence theorem, the integral of [tex]\vec F[/tex] across [tex]S[/tex] is equivalent to the integral of [tex]\mathrm{div}(\vec F)[/tex] over the interior of [tex]S[/tex], so that
[tex]\displaystyle \iint_S \vec F\cdot d\vec S = \iiint_{\mathrm{int}(S)} \mathrm{div}(\vec F)\,dV[/tex]
The paraboloid meets the [tex]x,y[/tex]-plane in a circle with radius 3, so we have
[tex]\mathrm{int}(S) = \left\{(x,y,z) \mid x^2+y^2\le3 \text{ and } 0 \le z \le 9-x^2-y^2\right\}[/tex]
and
[tex]\displaystyle \iiint_{\mathrm{int}(S)} \mathrm{div}(\vec F) \,dV = \int_{-3}^3 \int_{-\sqrt{9-x^2}}^{\sqrt{9-x^2}} \int_0^{9-x^2-y^2} 6(x^2+y^2)\,dz\,dy\,dx[/tex]
Convert to cylindrical coordinates, with
[tex]\begin{cases}x = r\cos(\theta) \\ y = r\sin(\theta) \\ z = \zeta \\ dV = dx\,dy\,dz = r\,dr\,d\theta\,d\zeta\end{cases}[/tex]
so that [tex]x^2+y^2=r^2[/tex], and the domain of integration is the set
[tex]\left\{(r,\theta,\zeta) \mid 0 \le r \le 3\text{ and } 0 \le \theta\le2\pi \text{ and } 0 \le \zeta \le 9-r^2\right\}[/tex]
Now compute the integral.
[tex]\displaystyle \int_0^3 \int_0^{2\pi} \int_0^{9-r^2} 6r^2\cdot r\,d\zeta\,d\theta\,dr = 12\pi \int_0^3 \int_0^{9-r^2} r^3\, d\zeta \, dr \\\\ ~~~~~~~~~~~~ = 12\pi \int_0^3 r^3 (9 - r^2) \, dr \\\\ ~~~~~~~~~~~~ = 12\pi \int_0^3 (9r^3 - r^5) \, dr \\\\ ~~~~~~~~~~~~ = 12\pi \left(\frac94 r^4 - \frac16 r^6\right)\bigg|_0^3 = 12\pi \left(\frac94\cdot3^4-\frac16\cdot3^6\right) = \boxed{729\pi}[/tex]