where $S$ describes the six faces of the cube. On these six faces, we have \begin{alignat*}{4} S_1 &:\quad \mb{F} &&= [0, 0, z^2 - 1] \qquad &\mb{n} &= [0, -1, 0] \\ S_2 &:\quad \mb{F} &&= [x^2(x-1), 0, z^2 - 1] \qquad &\mb{n} &= [0, -1, 0] \\ S_3 &:\quad \mb{F} &&= [0, 0, z^2 - 1] \qquad &\mb{n} &= [-1, 0, 0] \\ S_4 &:\quad \mb{F} &&= [0, y(y-1)^2, z^2 - 1] \qquad &\mb{n} &= [1, 0, 0] \\ S_5 &:\quad \mb{F} &&= [(x-1)x^2 y, (y-1)^2 xy, - 1] \qquad &\mb{n} &= [0, 0, -1] \\ S_6 &:\quad \mb{F} &&= [(x-1)x^2 y, (y-1)^2 xy, 0] \qquad &\mb{n} &= [0, 0, 1] \end{alignat*}

where you should be able to figure out which face corresponds to which $S_i$ by the normal. Thus, we have \begin{alignat*}{4} S_1 &:\quad \mb{F} \cdot\mb{n} &&= 0 \\ S_2 &:\quad \mb{F} \cdot\mb{n}&&= 0 \\ S_3 &:\quad \mb{F} \cdot\mb{n}&&= 0 \\ S_4 &:\quad \mb{F} \cdot\mb{n}&&= 0 \\ S_5 &:\quad \mb{F} \cdot\mb{n}&&= 1 \\ S_6 &:\quad \mb{F} \cdot\mb{n}&&= 0. \end{alignat*}

and so the right hand side of the divergence theorem should be \[ \iint_{S_5} (1) \, \de{S} = 1, \]

since the area of the face is 1. Next we compute the left hand side. This gives: \[ \int_0^1 \int_0^1 \int_0^1 [x + 2 + xy(x-2) + 3xy^2 + 2z] \, \de{x} \de{y} \de{z} = 1, \]

after some work.