At the bottom of the cylinder, $z=0$. Since $F_z(x,y,0)=0$ we find
$$\begin{align}
\int_{S_1} \vec F \cdot \hat n\, dS&=\int_{S_1}\vec F(x,y,0)\cdot \hat z\,dS\\\\
&=0
\end{align}$$
and there is no flux contributed.
At the top of the cylinder, $z=5$. Since $F_z(x,y,5)=15$ we find
$$\begin{align}
\int_{S_2} \vec F \cdot \hat n\, dS&=\int_{S_2}\vec F(x,y,5)\cdot \hat z\,dS\\\\
&=\int_0^{2\pi}\int_0^2 (15)\,\rho\,d\rho\,d\phi\\\\
&=60\pi
\end{align}$$
Putting it all together, we find
$$\begin{align}
\int_S \vec F\cdot \hat n \,dS&=\int_V \nabla \cdot \vec F\,dV-\int_{S_2} \vec F \cdot \hat n\, dS-\int_{S_1} \vec F \cdot \hat n\, dS\\\\
&=120\pi-60\pi-0\\\\
&=60\pi
\end{align}$$
First of all, it is always good practice to choose a parametrization that represents the whole surface.
Secondly, $-12 \pi$ is incorrect answer unless there are typo in the question.
Thirdly, there must be mistakes how you arrived at the final integral. It does not simplify to what you wrote. Also note that if you do the dot product with unit normal vector, you will have to use a factor to map surface area element to area of the projection in yz plane. Here is the correct integral using your approach,
The unit normal vector is $ \hat n = \left(\cfrac{x}{2}, \cfrac{y}{2}, 0 \right)$
So, $\vec F \cdot \hat n = x^2 (y + z) + \cfrac{y^3}{2}$
Now using $x = \sqrt{4-y^2}$ and $ds = \sqrt{1 + (x')^2} \ dy \ dz$, the integral is,
$\displaystyle 2 \int_{-2}^2 \int_{y-2}^0 \left((4-y^2) (y + z) + \frac{y^3}{2} \right) \frac{2}{\sqrt{4-y^2}} \ dz \ dy$
$ = - 40 \pi$
In fact, you can confirm this result using divergence theorem.
Here is what would be a better parametrization given it is a cylinder,
$r(\theta, z) = (2\cos\theta, 2 \sin \theta, z), 0 \leq \theta \leq 2\pi, 2 \sin\theta - 2 \leq z \leq 0$
$r_{\theta} \times r_z = (2 \cos\theta, 2\sin\theta, 0)$
$ \vec F(x,y,z)=(xy+2xz,x^2+y^2,xy-z)$
So, $\vec F (r(\theta, z)) = (2 \sin2\theta + 4z \cos\theta, 4, 2 \sin2\theta - z)$
$\vec F \cdot (r_{\theta} \times r_z) = 4 \sin2\theta \cos\theta + 8 z \cos^2\theta + 8 \sin\theta$
So integral is,
$\displaystyle \int_0^{2\pi} \int_{2 \sin\theta - 2}^0 (4 \sin2\theta \cos\theta + 8 z \cos^2\theta + 8 \sin\theta) \ dz \ d\theta = - 40 \pi$
After you integrate wrt $z$ and while integrating wrt $\theta$, you can avoid integrating functions whose integral is zero over $(0, 2\pi)$. That should simplify working to some extent.
Best Answer
Note that the closed surface $S$ consists out of $3$ pieces: the cylinder mantle, the bottom disc and the upper disc. So you need to calculate three integrals.
Cylinder mantle: I'll set up one integral as an example. A parametrisation is $$ \Sigma: \begin{cases} x = R \cos t, \\ y = R \sin t, \\ z = u \end{cases} \qquad \qquad \text{with $t\in[0,2\pi]$ and $u \in [0,H]$.} $$ The normal $n$ is equal to the cross product of the partial derivatives: $$ n(t,u) = \frac{\partial \Sigma}{\partial t} \times \frac{\partial \Sigma}{\partial u} = \begin{bmatrix} - R \sin t \\ R \cos t \\ 0 \end{bmatrix} \times \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} = \begin{bmatrix} R \cos t \\ R \sin t \\ 0 \end{bmatrix}. $$ This is an outward point normal vector. Please note that $n$ is not a unit vector! Then the inner product becomes $F(t,u)\cdot n(t,u)= R^2 \cos^2 t u + R^3 \cos t \sin^2 t$. The integral over the cylinder becomes $$ \iint_{\Sigma_1} F\cdot n = \int_0^{2\pi} \int_0^H R^2 \cos^2 t u + R^3 \cos t \sin^2 t \,du dt. $$
Upper and bottom disc: For the upper/bottom disc you need to find a parametrisation and set up an integral as hereabove. You need to make sure that the normal vector points upward on the upper disc, i.e., the $z$-coordinate needs to be positive. (For the bottom disc the $z$-coordinate of the normal vector needs to be negative.)
Can you take it from here?