a) Flux through $S_1$ (top disk in plane $y+z=2$):
Mistake: you multiplied by $\sqrt2$. Dot product with normal vector $(0, 1, 1)$ takes care of the factor for surface area element. If you are multiplying by $\sqrt2$, then you also need to make sure you have dot product with unit normal vector which is $\hat{n}=\frac{1}{\sqrt2}(0, 1, 1)$.
Parametrization of the surface: $(r \cos t, r \sin t, 2 - r \sin t), 0 \leq r \leq 1, 0 \leq t \leq 2\pi$
and $\vec{n} = (0, 1, 1)$.
$I_{S_1} = \displaystyle \int_0^{2\pi} \int_0^1 (4 r \cos t, -3r \sin t, 4 - 2 r \sin t) \cdot (0, 1, 1) \ r \ dr \ dt$
b) Flux through $S_2$ (cylindrical surface $x^2+y^2 = 1$):
Mistake: you hade a mistake in your normal vector and that is doubling the flux.
Parametrization of the surface: $s(z, t)=(\cos t, \sin t, z), 0 \leq t \leq 2\pi, 0 \leq z \leq 2 - \sin t$
and $\vec{n} = (x, y, 0) = (\cos t, \sin t, 0)$.
$I_{S_2} = \displaystyle \int_0^{2\pi} \int_0^{2-\sin t} (4 \cos t, -3 \sin t, 2z) \cdot (\cos t, \sin t, 0) \ dz \ dt$
[For outward normal vector to the cylinder it is easy to see that from axis of the cylinder $(0, 0, z)$ to a point on the cylinder orthogonal to the surface $(x, y, z)$, the vector is $(x, y, 0)$. Otherwise you can take partial derivatives of $s(z, t)$ and do the cross product $s'_t \times s'_z$ and that would give you the same.]
I've figured out my mistake, thanks to @Event Horizon. My impression from the Help Center page is that I shouldn't delete my question, so I'll outline what went wrong:
I should've applied the divergence theorem to $\nabla \times \mathbf{F}$ instead, so that the statement of the theorem becomes
$$ \iiint_\Omega \nabla \cdot (\nabla \times \mathbf{F}) = \iint_S \nabla \times \mathbf{F} \cdot d\mathbf{S} + \iint_\Sigma \nabla \times \mathbf{F}\cdot d\mathbf{S}. $$
The left is $0$ by since the divergence of a curl is $0$. From $\nabla \times \mathbf{F} = (z^2+x,0,-z-3)$, the correct computation is
$$ \iint_\Sigma \nabla \times \mathbf{F}\cdot d\mathbf{S} = \iint_{x^2+y^2\leq 4 \\ z=2} (-z-3 )dA = \iint-5dA = -20\pi$$ which matches what I got using Stokes' theorem (up to sign - but this just depends on the orientation of $S$).
Best Answer
Your computation is correct. Indeed, the direct computation gives $$ \begin{align}\int_S \mathbf{F\cdot n}\ dS&=\iint_{x^2+y^2\leq 4^2}(x,y,4^2)\cdot(0,0,1) dx dy \\&\;+\iint_{x^2+y^2\leq 4^2}(x,y,x^2+y^2)\cdot(\frac{x}{\sqrt{x^2+y^2}},\frac{y}{\sqrt{x^2+y^2}},-1) dx dy \\&=16\cdot 16\pi+\int_{\theta=0}^{2\pi}\int_{r=0}^4(r-r^2)r dr d\theta\\ &=256\pi +2\pi\left[\frac{r^3}{3}-\frac{r^4}{4}\right]_0^4 =256\pi-\frac{256\pi}{3}=\frac {512\pi}{3}. \end{align}$$
P.S. In the book's answer (see the comment below), the volume of the cone is $\frac{\pi}{3}\cdot 4^2\cdot 4$ which is equal to $\frac{64\pi}{3}$, it is not $\frac{16\pi}{3}$. Therefore a factor $4$ is missing, and $4\cdot \frac{ 128\pi}{3}$ gives you the correct answer $\frac{512\pi}{3} $.