[tex]\vec F(x,y,z)=(2x^3+y^3)\,\vec\imath+(y^3+z^3)\,\vec\jmath+3y^2z\,\vec k[/tex]
has divergence
[tex]\mathrm{div}\vec F(x,y,z)=6x^2+3y^2+3y^2=6(x^2+y^2)[/tex]
Then by the divergence theorem, the flux of [tex]\vec F[/tex] across [tex]S[/tex] is equal to the integral of [tex]\mathrm{div}\vec F[/tex] over the interior of [tex]S[/tex]. In cylindrical coordinates, this integral is
[tex]\displaystyle6\int_0^{2\pi}\int_0^1\int_0^{1-r^2}r^3\,\mathrm dz\,\mathrm dr\,\mathrm d\theta=\boxed{\pi}[/tex]
