Consider the vector field $\mathbf{F}(\mathbf{x})=\left(\left(3 x^{3}-x^{2}\right) y,\left(y^{3}-2 y^{2}+y\right) x, z^{2}-1\right)$ and let $S$ be the surface of a unit cube with one corner at $(0,0,0)$, another corner at $(1,1,1)$ and aligned with edges along the $x$-, $y$ - and $z$-axes. Use the divergence theorem to evaluate

$$I=\int_{S} \mathbf{F} \cdot d \mathbf{S}$$

Verify your result by calculating the integral directly.