Recall that by divercence theorem we obtain the flux through the whole surface enclosing the volume but in that case we are looking for the flux through the cylinder and bases surfaces included in the first octant.
Indeed by divergence theorem we obtain:
$$F=\int_V -3y^2 dV=\int_{0}^{\pi/2}\int_{0}^{a}\int_{0}^{b}-3r^3\sin^2\theta \:dz\:dr\:d\theta=-\frac{3a^4b\pi}{16}$$
And by direct calculation for the whole surface enclosing the volume:
$$F_1=\int_{S_1} \vec {F}\cdot (-\vec i)dS=0$$
$$F_{2,1}=\int_{S_{2,1}} \vec {F}\cdot \vec n \;dS
=\int_{S_{2,1}} (z\hat{i} +x\hat{j})\cdot \frac{x\hat{i} +y\hat{j}}{a}\;dS
=\int_{S_{2,1}} (xz+xy)\;dS=
\\=\int_{0}^{\pi/2}\int_{0}^{b} (az\cos \theta+a^2\cos\theta \sin \theta) \;dz\;d\theta =\int_{0}^{\pi/2} \left(\frac12ab^2\cos \theta+a^2b\cos\theta \sin \theta\right) \;d\theta =
\\=\frac12ab^2+\frac12a^2b$$
$$F_{2,2}=\int_{S_{2,2}} \vec {F}\cdot (-\vec {j}) \;dS
=\int_{S_{2,2}} -x \;dS=\int_{0}^{a}\int_{0}^{b}-x \;dz\;dx=-\frac12a^2b$$
$$F_{2,3}=\int_{S_{2,3}} \vec {F}\cdot (-\vec {i}) \;dS
=\int_{S_{2,2}} -z \;dS=\int_{0}^{a}\int_{0}^{b}-z \;dz\;dx=-\frac12ab^2$$
$$F_3=\int_{S_3} \vec {F}\cdot \vec k\;dS=\int_{S_3} -3y^2b\;dS=-3\int_{0}^{\pi/2}\int_{0}^{a}r^3b\sin^2 \theta\;dr\;d\theta =-\frac{3a^4b\pi}{16}$$
$$F=F_1+F_2+F_3=-\frac{3a^4b\pi}{16}$$
which agrees with the result obtained by divergenge theorem.
The divergence theorem is about a three-dimensional body $B$, its boundary surface $\partial B$, and a vector field ${\bf v}$ with domain $\Omega\supset B$. The theorem says that
$$\int_{\partial B}{\bf v}\cdot{\bf n}\>dS=\int_B{\rm div}\,{\bf v}({\bf x})\> dV\ .\tag{1}$$
In your case $B$ is a half ball, and ${\bf v}={\rm curl}(\vec F)$. It follows that ${\rm div}\,{\bf v}({\bf x})\equiv 0$ in $B$. But $\partial B$ is not just the spherical part of $\partial B$. The boundary of $B$ also contains a plane circular disc in the $(x,y)$-plane. It follows that the right hand side of $(1)$ consists of two terms, one of which you have already computed. Computing the flow of ${\bf v}$ through the disc is simpler, since ${\bf n}\equiv(0,0,-1)$ and $dS=dA$, where $dA$ is the standard area element in the $(x,y)$-plane.
Best Answer
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}$$