Here's the task I'm a bit confused with:
Find using Gauss's theorem the outward flux seen from the point $\left(\begin{matrix}0\\0\\1\end{matrix}\right)$ of the vector field $$\mathbf{v}=\left(\begin{matrix}x^2yz\\xy^2z\\z\end{matrix}\right)$$ through the area $M$ defined by $$M:=\left\{ (x,y,z)\in\mathbb{R}^3:x^2+y^2=z,\;0\leq z \leq 1\right\}$$
I don't understand the meaning of the "seen from the point". Does it influence somehow the solution?
Here's how I do this.
At first I find divergence of $\mathbf{v}$: $$\operatorname{div}(\mathbf{v})=2xyz+2xyz+1=4xyz+1$$ To find the volume integral I make the following transformations: $$-1\leq x\leq 1\\ -\sqrt{1-x^2}\leq y \leq \sqrt{1-x^2}\\0\leq z \leq 1$$
Applying divergence theorem we obtain $$\int_{\partial V}\mathbf{v}\cdot\mathbf{n}\operatorname{do}=\int_V\operatorname{div}(\mathbf{v})d\mu=\int_{-1}^1dx\int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}}dy\int_{0}^14xyz+1dz\\=\dots\\=\pi$$
I skipped the routine calculation of the volume integral which I've done on the paper.
Does it seem to be correct? Is there anything to be done with this additional "seen from the point" condition?
It seems that your boundaries of integration are incorrect. The domain ${M = \{(x,y,z) \,|\, x^{2} + y^{2} = z , 0 \leq z \leq 1 \}}$ corresponds qualitatively to the surface of a cone of height $1$ oriented along the $z-axis$, without its top-hat. Indeed, one should note that the top of this cone, which corresponds to the surface $\{(x,y,1) \,|\, x^{2} + y^{2} \leq 1 \}$ does not belong to the domain $M$.
If we want to determine the flux of the vector field $\boldsymbol{v}$ through the surface $M$ using Gauss's Theorem one may proceed as follows
$$ \begin{aligned} \int_{\rm M} \boldsymbol{v} \cdot \boldsymbol{n} \, d S & = \int_{\partial \rm Cone} \boldsymbol{v} \cdot \boldsymbol{n} \, d S - \int_{\rm Top} \boldsymbol{v} \cdot \boldsymbol{n} \, d S \\ & = \int_{\rm Cone} \text{div} (\boldsymbol{v}) \, d V - \int_{\rm Top} \boldsymbol{v} \cdot \boldsymbol{n} \, d S \, , \end{aligned}$$
where $\int_{M} d S$ corresponds to the outgoing flux through the surface $M$, $\int_{\partial \rm Cone} d S $ corresponds to the flux through the surface of the entire cone with its top, and $\int_{\rm Top} d S$ corresponds to the flux through the top of the cone.
The integration $\int_{\partial \rm Cone} d S$ can be estimated thanks to Gauss's Theorem as follows. The domain of integration is given by
$$\begin{cases} 0 \leq z \leq 1 \, , \\ -\sqrt{z} \leq x \leq \sqrt{z} \, , \\ - \sqrt{z - x^{2}} \leq y \leq \sqrt{z - x^{2}} \, . \end{cases} $$
(The boundaries are different than yours, since you are integrating on the entire cylinder) so that it becomes
$$ \begin{aligned} \int_{\rm Cone} \text{div} (\boldsymbol{v}) \, d V & = \int_{0}^{1} d z \int_{-\sqrt{z}}^{\sqrt{z}} dx \, \int_{-\sqrt{z - x^{2}}}^{\sqrt{z - x^{2}}} d y \, (4 x y z + 1) \\ & = \frac{\pi}{2} \end{aligned}$$
On the other hand, the flux through the top $\int_{\rm Top} d S$ is given by
$$ \begin{aligned} \int_{\rm Top} \boldsymbol{v} \cdot \boldsymbol{n} \, d S & = \int_{-1}^{1} d x \, \int_{- \sqrt{1 - x^{2}}}^{\sqrt{1-x^{2}}} d y \, \, (x^{2} y \,,\, x y^{2} \,,\, 1) \cdot (0,0,1) \\ & = \int_{-1}^{1} d x \, \int_{- \sqrt{1 - x^{2}}}^{\sqrt{1-x^{2}}} d y \, 1 \\ & = \pi \end{aligned} $$ (Note the fact that on the top, the outgoing normal is given by $(0,0,+1)$)
As a conclusion, we obtain that $\int_{M} \boldsymbol{v} \cdot \boldsymbol{n} \, dS = - \frac{\pi }{2}$