Gauss's divergence integration

I'm not sure on how to compute this integral using Gauss's divergence theorem. Can someone please explain.

$\iint_S \vec{F} \cdot \vec{n}\quad dS$

$\vec{n}$ is outward normal and S is exterior surface of the cylinder $x^2 + y^2 ≤ 1, 0 ≤ z ≤ 1 \quad$ and $\quad \vec{F} (x, y, z) = (x^2y, z − xy^2, z^2)$

$\quad \vec{F} (x, y, z) =(x^2y, z − xy^2, z^2)$

$\nabla·\vec F=2xy-2xy+2z=2z$

$$\int_S\vec{F}.\hat nd\mathbf s=\int_V\nabla·\vec F\mathrm dv=\int_0^1\int_{-1}^1\int_{-\sqrt{1-y^2}}^{\sqrt{1-y^2}}2z\mathrm dx\mathrm dy\mathrm dz=$$

$$=\int_0^1\int_{-1}^1\left[2zx\right]_{x=-\sqrt{1-y^2}}^{x=\sqrt{1-y^2}}\mathrm dy\mathrm dz=\int_0^1\int_{-1}^12z\left(\sqrt{1-y^2}-\left(-\sqrt{1-y^2}\right)\right)\mathrm dy\mathrm dz$$

$$=\int_0^1\int_{-1}^1\left(4z\sqrt{1-y^2}\right)\mathrm dy\mathrm dz=\int_0^14z\left[y\sqrt{1-y^2}+\arcsin y\right]_{-1}^1\mathrm dz=$$

$$=\int_0^12\pi z\mathrm dz=2\pi[z^2/2]_0^1=\pi$$

The theorem you wrote is incorrect. it should be $$\iiint_V\nabla .\vec{F}dV=\iint_S\vec{F}.\hat ndS$$ where i suggest you evaluate the lhs so you do not need to calculate $\hat n$