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.$$
Based on the fact that $z = \sqrt{4-y^2}$ :
$$
\sigma(r,\theta) = (x,y,z) = (r\cos\theta, r\sin\theta, \sqrt{4-(r\sin\theta)^2}) \\
r \in [0,1],\ \theta \in [0,2\pi[
$$
According to WA, (And my understanding of the surface), this seems right.
The cross product is painful, I agree...
$$
\sigma_r = (\cos\theta, \sin\theta, -\frac{r\sin^2\theta}{\sqrt{4-r^2\sin^2\theta}}) \\
\sigma_\theta = (-r\sin\theta, r\cos\theta, -\frac{r^2\sin\theta\cos\theta}{\sqrt{4-r^2\sin^2\theta}}) \\
\sigma_r \times \sigma_\theta = (0,\frac{r^2(\cos^2\theta\sin\theta + \sin^3\theta)}{\sqrt{4-r^2\sin^2\theta}},r)
$$
Regarding your parametrization it doesn't work for $\theta < \pi/3, \theta > 2\pi/3$, the square roots in the bounds give complex numbers but I see what you meant to do. The problem is that the bound for $x$ doesn't define a circle of radius $1$ as it should be. Here you are mixing the cylinder of radius $2$ (which gives you the $y=2\cos\theta$) and the equations $x = \pm \sqrt{1-y^2}$ and to make this equation work $y$ needs to be between $-1$ and $1$ but by its nature its range is $[-2,2]$.
Checking Stokes theorem with this surface:
$$
A= \iint_S rot \vec{F} \cdot d\vec{n} = \iint_S (0,x,1) \cdot d\vec{n} =
\iint_S\left( r\cos\theta \cdot \frac{r^2(\cos^2\theta\sin\theta + \sin^3\theta)}{\sqrt{4-r^2\sin^2\theta}} + r \right)dS = \\
\underbrace{\int_0^1 dr \int_0^{2\pi} d\theta \frac{r^3(\cos^3\theta\sin\theta + \sin^3\theta\cos\theta)}{\sqrt{4-r^2\sin^2\theta}}}_{=I_1} + \int_0^1 dr \int_0^{2\pi} d\theta\ r = 0 + \pi
$$
Then
$$
\partial S(\theta) = (\cos\theta,\sin\theta, \sqrt{4-\sin^2\theta}), \ \partial S'(\theta) = (-\sin\theta,\cos\theta, -\frac{\sin\theta \cos\theta}{\sqrt{4-\sin^2\theta}}),\ \theta\in [0,2\pi[ \\
B = \oint_{\partial S} \vec{F}\cdot d\vec{l} = \oint_{\partial S} (zx-y,0,0) \cdot d \vec{l} = \int_0^{2\pi} (\cos\theta \sqrt{4-\sin^2\theta} -\sin\theta) (-\sin\theta)\; d\theta = \pi
$$
Therefore $A = B$. However I agree that the integrals are not be that simple but clearly doable which may indicate that the parametrization I chose was maybe not the best, but at least it worked !
$$
I_1 = \int_0^1 dr\ r^3 \int_0^{2\pi}d\theta \cos\theta \frac{(1-\sin^2\theta)\sin\theta + \sin^3\theta}{\sqrt{4-r^2\sin^2\theta}} \overset{w = \sin\theta, dw = \cos\theta d\theta}{=} \int_0^1 dr\ r^3\int_0^0 dw \frac{w}{\sqrt{4-r^2 w^2}} = 0
$$
Best Answer
You just need to solve your inequality for $r$: \begin{align*} a^2 &\geq 2 r^2 \cos^2(\theta) + 5 r^2 \sin^2(\theta) + 4 r^2 \sin(\theta) \cos(\theta)\\ &= r^2(2 \cos^2(\theta) + 5 \sin^2(\theta) + 4 \sin(\theta) \cos(\theta))\\ &= r^2(2 + 3 \sin^2(\theta) + 4 \sin(\theta) \cos(\theta)) \, . \end{align*} Then $$ r^2 \leq \frac{a^2}{2 + 3 \sin^2(\theta) + 4 \sin(\theta) \cos(\theta)} $$ (Note that the denominator is always $\geq 0$, so the inequality is preserved.) Thus your bounds of integration should be \begin{align*} &0 \leq \theta \leq 2 \pi\\ &0 \leq r \leq \frac{a}{\sqrt{2 + 3 \sin^2(\theta) + 4 \sin(\theta) \cos(\theta)}} \, . \end{align*}