Rigorous proof of Gauss divergence theorem 
Gauss divergence theorem: If $V$ is a compact volume, $S$ its boundary being piecewise smooth and $F$ is a continuously differentiable vector field defined on a neighborhood of $V$, then we have:
$$\iiint_V \left(\mathbf{\nabla}\cdot\mathbf{F}\right)\,dV=
\unicode{x222F} (\mathbf{F}\cdot\mathbf{n})\,dS $$

Right now I am taking a real analysis course. The lecturer discusses the proof of Stokes curl theorem but not Gauss divergence theorem. Can someone provide a rigorous proof of Gauss divergence theorem so that a real analysis student can understand?

If writing the proof is tedious, please provide a link of rigorous proof of Gauss divergence theorem

 A: Theorem:

Let R be the closure of a bounded open set in $\Bbb R^3$ with its boundary $\partial R$ being a simple surface. Suppose $F(x,y,z)=(a(x,y,z),b(x,y,z),c(x,y,z)),$ where the functions are continuous in R and have continuous and bounded first derivatives in the interior points of R. Then $\iiint_R div(F)dxdydz=\iint_{\partial R} F\cdot n \;dA.$

Proof:
First, we shall prove the following.

If every point $Q\in R$ has a neighborhood $B_Q(\epsilon_Q)$ such that the formula holds for all $a,b,c$ having support in $B_Q(\epsilon_Q)$, it holds for general $a,b,c$.

Define $\psi_Q(P)=\begin{cases}(\epsilon_Q^2-4\overline {PQ}^2)^2\;\text{for}\;\overline {PQ}\lt\frac12\epsilon_Q\\ 0\;\text{for}\;\overline {PQ}\ge\frac12\epsilon_Q\end{cases}$ which are continuous and have continuous first derivatives for all P. Since R is compact, by Heine-Borel it has a finite subcover $\{B_{Q_1}(\frac 12\epsilon_{Q_1}),...,B_{Q_N}(\frac 12\epsilon_{Q_N})\}$. We then introduce functions $\chi_i(P)=\dfrac{\psi_{Q_i}(P)}{\psi_{Q_1}(P)+...+\psi_{Q_N}(P)}\in C^1(R)$  satisfying the conditions for a partition of unity, i.e. $\chi_i(P)\ge 0, \sum_i \chi_i(P)=1,$ and $\chi_i(P)$ have support in $B_{Q_i}(\frac 12\epsilon_{Q_i})$. We may then write $a=\sum_i a\chi_i(P)$, where the individual terms are supported in $\{B_{Q_i}(\frac 12\epsilon_{Q_i})$. Decomposing $b,c$ similarly, we see that since the formula holds for individual terms, it holds for the whole expression.
Hence, we only have to prove the theorem for $a,b,c$ supported in an arbitrarily small neighborhood of a point Q.
For a point Q interior to R, $a=b=c=0$ on the boundary, while if we pick a cube $I\supset R$, we have $\iiint_R c_z dxdydz= \iiint_I c_z dxdydz$, which vanishes since the boundaries of integration are outside of the support of $c$. Similar arguments may be made for $a_x,b_y$, and thus the left hand side of the formula is also 0, which proves the theorem for this case.
For $Q\in \partial R$, we may choose a suitable $\epsilon_Q$ such that the points near Q may be parametrized by $\Bbb X: U(u,v)\rightarrow B_Q(\epsilon_Q)\cap \partial R$. By a suitable rigid motion of the space (this doesn't alter the formula, as is easily verified), we may assume none of the components of $\Bbb X_u \times \Bbb X_v$ vanishes at Q. Thus, for a refined and perhaps smaller $\epsilon_Q$, the component does not vanish on $B_Q(\epsilon_Q)$, and one may express any one of the three vairables $x,y,z$ as a function of the other two. For example, we may write $z=\phi(x,y)$ for points sufficiently near Q.
Consider $(x,y,z)$ near the point Q. There are two cases: $(1):z\gt \phi(x,y)$ or $(2):z\lt \phi(x,y)$. If (1) is inside R  while (2) is outside R, we may pick a parametrization $\Bbb X=(v,u,\phi(u,v))$. If (1) is outside R  while (2) is inside R, we may pick instead $\Bbb X=(u,v,\phi(u,v))$.
WLOG, assume the second case. Then the normal vector $n=\dfrac{\Bbb X_u \times \Bbb X_v}{|\Bbb X_u \times \Bbb X_v|}=(\dfrac{-\phi_u}{\sqrt{1+\phi_u^2+\phi_v^2}},\dfrac{-\phi_v}{\sqrt{1+\phi_u^2+\phi_v^2}},\dfrac{1}{\sqrt{1+\phi_u^2+\phi_v^2}})$. We denote these three components by $\xi,\eta,\zeta$ respectively, and consider the last component.
Since the area element $dA=\sqrt{1+\phi_u^2+\phi_v^2}, \iint_{\partial R} c\zeta dA= \iint_U c dudv$. By choosing the orientation $x=u,y=v,z=\phi(u,v)$, we see that $\iint_U c dudv=\iint_U c dxdy$. On the other hand, we may extend c continuously as 0 outside of its domain. Then, $$\iiint_R c_z dxdydz=\iiint_{z\le \phi(x,y)} c_z dxdydz=\iint(\int_{-\infty}^{\phi(x,y)} c_z dz)dxdy=\iint (c(x,y,\phi(x,y))-0) dxdy.$$
Therefore, the two integrals have the same value. Again, similar argument may be made for the other two variables, which completes the proof.
My professor demonstrated this proof in the Calculus II course I'm currently taking. Since you've stated you're a real analysis student, I skipped over some definitions as you have probably encountered them before. Hope this version of divergence theorem is sufficient for your application needs.
