Let's call $S_c$ the curved surface of the cylindrical solid of radius $2$ about the interval $[0,3]$ along the $z$-axis, $S_t$ the top flat surface, $S_b$ the bottom flat surface, and similarly index the respective outward unit normal vectors. We are trying, then, to find $\int_{S_c}\mathbf{F}\cdot\mathbf{n}_c\,dS_c.$
Note that we have $\mathbf{n}_c=\frac{\nabla(x^2+y^2-4)}{\lVert\nabla(x^2+y^2-4)\rVert}=\frac{2x\hat i+2y\hat j+0\hat k}{\lVert 2x\hat i+2y\hat j+0\hat k\rVert}=\frac{2}{\sqrt{(2x)^2+(2y)^2+0^2}}\left(x\hat i+y\hat j+0\hat k\right)=\frac{1}{\sqrt{x^2+y^2}}\left(x\hat i+y\hat j+0\hat k\right),$ and recalling that $x^2+y^2=4$ on $S_c$, we have $\mathbf{n}_c=\frac{1}{2}\left(x\hat i+y\hat j+0\hat k\right).$ Thus, $\mathbf{F}\cdot\mathbf{n}_c=\frac{1}{2}(x-1)y^2.$
Also, $\mathbf{n}_t=\hat k$, and $\mathbf{n}_b=-\hat k$, so $\mathbf{F}\cdot\mathbf{n}_t=(y^2\hat i-y\hat j+3xy\hat k)\cdot\hat k=3xy,$ and $\mathbf{F}\cdot\mathbf{n}_b=(y^2\hat i-y\hat j+0\hat k)\cdot\hat k=0.$
By divergence theorem, $\int_V(\nabla\cdot\mathbf{F})\,dV=\int_S(\mathbf{F}\cdot\mathbf{n})\,dS,$ where $V$ is the cylindrical solid described at the start and $S$ is the full surface of said solid. It is clear that $\int_S(\mathbf{F}\cdot\mathbf{n})\,dS=\int_{S_c}(\mathbf{F}\cdot\mathbf{n}_c)\,dS_c+\int_{S_t}(\mathbf{F}\cdot\mathbf{n}_t)\,dS_t+\int_{S_b}(\mathbf{F}\cdot\mathbf{n}_b)\,dS_b,$ so by the work above, $\int_S(\mathbf{F}\cdot\mathbf{n})\,dS=\int_{S_c}(\mathbf{F}\cdot\mathbf{n}_c)\,dS_c+3\int_{S_t}xy\,dS_t.$ Thus, since $\nabla\cdot\mathbf{F}=\left(\frac{\partial}{\partial x}\hat i+\frac{\partial}{\partial y}\hat j+\frac{\partial}{\partial z}\hat k\right)\cdot(y^2\hat i-y\hat j+xyz\hat k)=-x+xy=x(y-1),$ then we have $\int_{S_c}(\mathbf{F}\cdot\mathbf{n}_c)\,dS_c=\int_Vx(y-1)\,dV-3\int_{S_t}xy\,dS_t.\qquad (\#)$
Now, since we're dealing with a cylinder, it makes sense to use cylindrical coordinates to evaluate the volume integral. In particular, $dV=r\,dz\,dr\,d\theta$, $x=r\cos\theta$, $y=r\sin\theta$, $0\leq z\leq 3$, $0\leq r\leq 2$, and $0\leq\theta\leq 2\pi$--hopefully, that should be enough to let you evaluate $\int_Vx(y-1)\,dV$ as a triple integral.
We can readily parameterize $S_t$ by $\mathbf{r}(x,y)=x\hat i+y\hat j+3\hat k$, where $(x,y)$ varies over the region $C=\{(x,y)\in\Bbb R^2:x^2+y^2\leq 4\}$ in the $xy$-plane. Now, $dS_t=\left\lVert\frac{\partial\mathbf{r}}{\partial x}\times\frac{\partial\mathbf{r}}{\partial y}\right\rVert\,dA,$ so since $\frac{\partial\mathbf{r}}{\partial x}\times\frac{\partial\mathbf{r}}{\partial y}=\left|\begin{array}{ccc}\hat i & \hat j & \hat k\\ 1 & 0 & 0\\ 0 & 1 & 0\end{array}\right|=\hat k,$ so $dS_t=dA$ and $\int_{S_t}xy\,dS_t=\int_Cxy\,dA.$ As we are dealing with a circular region $C$, then polar coordinates will be the best way, with $dA=r\,dr\,d\theta$, $x=r\cos\theta$, $y=r\sin\theta$, $0\leq r\leq 2$, and $0\leq\theta\leq 2\pi$--again, hopefully this will be enough for you to calculate $\int_{S_t}xy\,dS_t$ as a double integral.
Once you've calculated the aforementioned triple and double integrals, the desired surface integral will follow from the equation $(\#)$.