Step 1: Skip the divergence theorem and add up the flux through each pair of opposite faces directly.
For \( \vec{F} = 4xz\,\hat{i} - y^2\,\hat{j} + yz\,\hat{k} \) over the unit cube, look at the \( x=0 \) and \( x=1 \) faces first. Only \( F_x = 4xz \) contributes here, and it is zero at \( x=0 \), so the net outward flux from this pair is:
\[ \int_0^1\!\!\int_0^1 4(1)z \, dy\, dz = 4 \times \tfrac{1}{2} = 2 \]
Step 2: Do the same for the \( y=0 \) and \( y=1 \) faces.
Only \( F_y = -y^2 \) contributes, giving:
\[ \int_0^1\!\!\int_0^1 \left[-(1)^2 - (-(0)^2)\right] dx\, dz = -1 \]
Step 3: Do the same for the \( z=0 \) and \( z=1 \) faces.
Only \( F_z = yz \) contributes, giving:
\[ \int_0^1\!\!\int_0^1 \left[y(1) - y(0)\right] dx\, dy = \int_0^1 y\, dy = \tfrac{1}{2} \]
Step 4: Final Answer.
Adding all three pairs:
\[ 2 - 1 + \tfrac{1}{2} = \boxed{\tfrac{3}{2}} \]
which matches the divergence theorem result. Therefore, option (D) is correct.