Let $D$ be the domain in space bounded by the planes $z = 0$ and $z = 2x$, along with the cylinder $x = 1 − y^2$, as graphed in the attached Figure, let $S$ be the boundary of $D$, and let $F = (x + y, y^2, 2z⟩$.
Verify the Divergence Theorem by finding the total outward flux of $F$ across $S$, and show this is equal to $$\int_V \text{div}(F) dV$$.
The surface S is piecewise smooth, comprising surfaces
- $S_1$ is part of the plane $z = 2x$,
- $S_2$, is part of the cylinder $x = 1 − y^2$ and
- $S_3$, is part of the plane $z = 0$.
I am unable to understand the following parametrization of the surfaces $S_1$, $S_2$ and $S_3$

Edit: My parametrization of $S_1$:
$$x=1-y^2, y=y, z=2(1-y^2)$$ with $y\in [-1,1]$, $x\in[0,1]$. Don't know it is correct or not?

No, your parametrization $S_1$ is not correct: $x=1-y^2$, $y=y$, $z=2(1-y^2)$ defines a curve for $y\in [-1,1]$ not a surface!
You may consider $S_1$ as a cartesian surface: $\mathbf{r}_1(x,y)=(x,y,2x)$ with $D_1=\{(x,y): x\in [0,1-y^2], y\in [-1,1]\}$. Therefore the outward flux through $S_1$ is $$\iint_{S_1}\mathbf{F}\cdot d\mathbf{S}=\iint_{D_1} \langle(x + y, y^2, 2(2x)),(-2,0,1)\rangle dx dy=\frac{16}{15}.$$
Similarly for $S_3$ let $\mathbf{r}_3(x,y)=(x,y,0)$ with $D_3=\{(x,y): x\in [0,1-y^2], y\in [-1,1]\}$ and we find $$\iint_{S_3}\mathbf{F}\cdot d\mathbf{S}=\iint_{D_3} \langle(x + y, y^2, 2(0)),(0,0,-1)\rangle dx dy=0.$$
As regards $S_2$, we may use $\mathbf{r}_2(y,z)=(1-y^2,y,z)$ with $D_2=\{(y,z): y\in [-1,1], z\in [0,2(1-y^2)]\}$ and therefore $$\iint_{S_2}\mathbf{F}\cdot d\mathbf{S}=\iint_{D_2} \langle((1-y^2) + y, y^2, 2z),(1,2y,0)\rangle dy dz=\frac{32}{15}.$$
Hence the total outward flux is $\frac{16}{15}+0+\frac{32}{15}=\frac{16}{5}$.
It remains to verify that the same result can be obtained as $$\iiint_D\text{div}(\mathbf{F})dV= \int_{y=-1}^1\int_{x=0}^{1-y^2}\int_{z=0}^{2x}(3+2y)dzdxdy.$$ Could you check this last step?