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.$$
Best Answer
The surface integral is better written as
$$\iint_S dS \, \hat{\mathbf{n}} \cdot \mathbf{F}$$
where $\hat{\mathbf{n}}$ is the unit surface normal. For surfaces of the form $z=f(x,y)$, the unit normal is
$$\hat{\mathbf{n}} = \frac{\displaystyle\left(-\frac{\partial f}{\partial x},-\frac{\partial f}{\partial y},1\right)}{\displaystyle\sqrt{\left(\frac{\partial f}{\partial x}\right)^2+\left(\frac{\partial f}{\partial y}\right)^2+1}}$$
Note that the normalizing factor in the denominator is exactly the stretching factor relating the surface area element to the differential $dx\,dy$:
$$dS = \sqrt{\left(\frac{\partial f}{\partial x}\right)^2+\left(\frac{\partial f}{\partial y}\right)^2+1} \, dx\, dy$$
Now, $\partial f/\partial x = 2$ and $\partial f/\partial y = -1$, so that
$$\iint_S dS \, \hat{\mathbf{n}} \cdot \mathbf{F} = \int_0^2 dx \, \int_0^2 dy \, [-(2) \cdot (3 x^2) - (-1) \cdot (-2 x y) + (1)\cdot (8)] $$
which is the result you wanted to see.