Compute the divergence:
[tex]\nabla\cdot\mathbf f(x,y,z)=\dfrac{\partial(x^2)}{\partial x}+\dfrac{\partial(xy)}{\partial y}+\dfrac{\partial z}{\partial z}=2x+x+1=3x+1[/tex]
By the divergence theorem, the flux is of [tex]\mathbf f[/tex] across [tex]\partial\mathcal E[/tex] (the boundary of the region [tex]\mathcal E[/tex]) is
[tex]\displaystyle\iint_{\partial\mathcal E}\mathbf f(x,y,z)\cdot\mathrm d\mathbf S=\iiint_{\mathcal E}\nabla\cdot\mathbf f(x,y,z)\,\mathrm dV[/tex]
We set up and compute the volume integral with respect to cylindrical coordinates.
[tex]\displaystyle\iiint_{\mathcal E}(3x+1)\,\mathrm dV=\int_{\theta=0}^{\theta=2\pi}\int_{r=0}^{r=2}\int_{z=0}^{z=4-r^2}(3r\cos\theta+1)r\,\mathrm dz\,\mathrm dr\,\mathrm d\theta=\frac{7\pi}2[/tex]