I haven't done a surface integral in a while so I am asking to get this checked.
$\mathbf{F} = \langle x, y, z\rangle$ and the surface is $z = xy + 1$ where $0\leq x,y\leq 1$.
$\hat{\mathbf{n}} = \nabla f/ \lvert\nabla f\rvert = \frac{1}{\sqrt{3}}\langle 1, 1, 1\rangle$
$dS = \frac{\lvert\nabla f\rvert dxdy}{\frac{\partial f}{\partial z}} = \sqrt{3}dxdy$
$\mathbf{F}\cdot\hat{\mathbf{n}} = \frac{1}{\sqrt{3}}(x+y+z) = \frac{1}{\sqrt{3}}(x+y+xy + 1)$
$$ \int_0^1\int_0^1(x + y + xy + 1)dxdy = \frac{9}{4} $$
However, when I used the divergence theorem, I obtained:
$$ \int_S(\mathbf{F}\cdot\hat{\mathbf{n}})dS = \int_V(\nabla\cdot\mathbf{F})dV $$ and $\nabla\cdot\mathbf{F} = 3$ so $$ \int 3dV = 3V = 3\frac{5}{4} = \frac{15}{4} $$
Which one is wrong or are both incorrect? If so, what is wrong?
I'm afraid neither is correct. You confused derivatives of $\mathbf F$ with calculating the surface normal. The Divergence Theorem does not apply, as you do not have a closed surface (so there's no volume $V$ enclosed). And to compute the surface integral you need to compute $$\int_0^1\int_0^1 \mathbf F(\mathbf g(x,y))\cdot \left(\frac{\partial\mathbf g}{\partial x}\times\frac{\partial\mathbf g}{\partial y}\right) \,dx\,dy\,,$$ where $\mathbf g(x,y) = (x,y,xy+1)$. This becomes $$\int_0^1\int_0^1 (1-xy)\,dx\,dy = \frac34\,.$$