I want to compute the flux of $F(x,y,z)=(xy)\pmb{i}+(yz)\pmb{j}+(zx)\pmb{k}$ across the surface $S$ given by the paraboloid $z=9-x^2-y^2$ over $[-1,1]\times[-1,1]$.
Ok, so I will compute $flux(F)=\int\int_S F\cdot N\;dS$.
Taking $f(x,y)=(x,y,9-x^2-y^2)$, then $$N={1\over\sqrt{1+4x^2+4y^2}} \begin{bmatrix} 2x \\ 2y \\ 1\end{bmatrix},$$ so $$N\cdot F dS={1\over\sqrt{1+4x^2+4y^2}} \begin{bmatrix} 2x \\ 2y \\ 1\end{bmatrix}\cdot\begin{bmatrix} xy \\ yz \\ zx\end{bmatrix}\sqrt{1+4x^2+4y^2}\;dA.$$
Thus, \begin{align} flux(F)&=\int\int_S F\cdot N\;dS=\int\int 2x^2y+2y^2z+zx\;dA\\ &=\int_{-1}^{1}\int_{-1}^{1}2x^2y+2y^2(9-x^2-y^2)+(9-x^2-y^2)x\;dy\;dx\\ &={968\over 45}. \end{align}
Right, but there is another to compute this using the divergence theorem. This is $$flux(F)=\int\int_S F\cdot N\;dS=\int\int\int_E\;div(F)\;dV.$$ We have that $div(F)=x+y+z$, so that \begin{align} flux(F)&=\int\int\int_E x+y+z\;dV=\\ &=\int_{-1}^{1}\int_{-1}^{1}x(9-x^2-y^2)+y(9-x^2-y^2)+{1\over 2}(9-x^2-y^2)^2\;dx\;dy\\ &={6266\over 45}. \end{align}
As you see, the results obtained with each method don't agree. I double-checked the computations by hand and with Mathematica, so I assume there must be a conceptual mistake or something like this.
Any help is appreciated.