The surface $S$ is the part of the plane $z=4+x+y$ that lies inside the cylinder $x^2+y^2=4$. We want to describe $S$ as a parametric surface $\vec r=\vec r(u,v)$. This is most easily accomplished in cylindrical coordinates, Where the equation of the plane takes the form $z=4+\rho\cos{\phi}+\rho\sin{\phi}=4+\rho(\cos{\phi}+\sin{\phi})$, and the equation of the cylinder takes the form $\rho=2$.
We have the following parametric representation of $S$:
$$\vec{r}(\rho,\phi)=\langle\rho\cos{\phi},\rho\sin{\phi},4+\rho(\cos{\phi}+\sin{\phi})\rangle,~~~0\leq\rho\leq2,0\leq\phi\leq2\pi.$$
Computing the partial derivatives and their cross product, we get
$$\vec{r}_\rho=\langle\cos{\phi},\sin{\phi},\cos{\phi}+\sin{\phi}\rangle$$
$$\vec{r}_\phi=\langle-\rho\sin{\phi},\rho\cos{\phi},\rho(\cos{\phi}-\sin{\phi})\rangle$$
$$\vec{r}_\rho\times\vec{r}_\phi=\langle-\rho,-\rho,\rho\rangle.$$
Thus, $\|\vec{r}_\rho\times\vec{r}_\phi\|=\sqrt{3}\rho$.
Rewriting the integrand in cylindrical coordinates, we have
$$f(x,y,z)=x^2z+y^2z=(x^2+y^2)z=\rho^2z.$$
Recall that on the surface $S$, we have $z=4+\rho(\cos{\phi}+\sin{\phi})$, so then
$$f(\vec{r}(\rho,\phi))=\rho^2(4+\rho(\cos{\phi}+\sin{\phi})).$$
Putting it all together, the surface integral evaluates to
$$I=\iint_{S}f(\vec{r})\,dS=\iint_{D}f(\vec{r}(\rho,\phi))\|\vec{r}_\rho\times\vec{r}_\phi\|\,dA\\
=\sqrt{3}\int_{0}^{2\pi}\int_{0}^{2}\rho^3(4+\rho(\cos{\phi}+\sin{\phi}))\,d\rho\,d\phi$$
We can make the evaluation of the above integral even easier by switching the order of integration and noting that $\cos\phi+\sin\phi$ is periodic in $\phi$ with period $2\pi$ and that the integral over one period is automatically zero. Hence, the value of the surface integral is:
$$I=\sqrt{3}\int_{0}^{2\pi}\int_{0}^{2}\rho^3(4+\rho(\cos{\phi}+\sin{\phi}))\,d\rho\,d\phi=\sqrt{3}\int_{0}^{2}\int_{0}^{2\pi}(4\rho^3+\rho^4(\cos{\phi}+\sin{\phi}))\,d\phi\,d\rho\\
=\sqrt{3}\int_{0}^{2}\int_{0}^{2\pi}4\rho^3\,d\phi\,d\rho\\
=2\sqrt{3}\pi\int_{0}^{2}4\rho^3\,d\rho\\
=2\sqrt{3}\pi(2^4)\\
=2^5\sqrt{3}\pi.$$
To make good use of the divergence theorem:
- verify that the surface integrals through the parts of the coordinate planes that, together with the two given surfaces, bound the region in the first octant, don't contribute to the flux (their share is $0$, which is easy to compute or even simply reason);
- note that due to symmetry in the two surfaces ($y \leftrightarrow z$) that make up $\partial T$ and with the given vector field, you expect the same flux through these two parts of $\partial T$.
By the divergence theorem you find a total flux of $\tfrac{8}{5}a^5$ so the flux through the part $S$ is half of this: $\tfrac{4}{5}a^5$.
You can verify this by actually computing the surface integral. To do this, we parametrize $S$ (the paraboloid $z=a^2-x^2$ cut off by the paraboloid $y=a^2-x^2$) as follows:
$$\vec r(u,v)=\left(u,v\left(a^2-u^2\right),a^2-u^2\right) \;,\; 0 \le u \le a\;,\; 0 \le v \le 1$$
Then:
$$\vec N=\frac{\partial \vec r}{\partial u}\times\frac{\partial \vec r}{\partial v}=\left(2a^2u-2u^3,0,4-u^2\right)$$
and:
$$\vec F \cdot \vec N=a^4-u^4$$
so:
$$\iint_S \vec F \cdot \vec N\,\mbox{d}S = \int_0^1 \int_0^a \left(a^4-u^4\right)\,\mbox{d}u \,\mbox{d}v=\frac{4}{5}a^5$$
Best Answer
Switching to a volume integral over the divergence indeed seems like a good idea, since now our problem is to simply find the volume of the region of integration (and then multiply it by $3$). Note that for a fixed value of $z$, say $z=0$, we get a 'slice' of the surface given by the equation $x^2+y^2=13$, which you may recognize as a circle with radius $\sqrt{13}.$ As $z$ increases, the circles become smaller, and when $z=13$ the slices disappear since the equation $x^2+y^2=13-z$ has no solutions for $z \geq 13$. So the lower $z$-limit of integration should be $4$ (which is given) and the upper limit should be $13$. Since the whole problem is cylindrically symmetric, cylindrical coordinates seem like a good idea. With this, the integral becomes $$3\int_4^{13} \int_0^{2\pi} \int_0^{\sqrt{13-z}} r\,dr\,d\theta\,dz,$$ since for any fixed $z$ the radius of the circle is $\sqrt{13-z}.$ This integral is not too hard to evaluate.
EDIT: As Ted pointed out, the above would be valid if the bottom was part of the surface integral, but it isn't. The fix to this is to subtract the flux integral over the bottom part, which is $$\int_{x^2+y^2=9,\\z=4} (x,y,z) \cdot (0,0,-1) \,dx \,dy = -36\pi.$$