Let $f: \mathbb{R}^n \supset \Omega \to \mathbb{R}$ and $u:\Omega \rightarrow \mathbb{R}$. For $$\mathrm{div} \Bigg (\frac{\nabla u}{\sqrt{1+|\nabla u|^2}}\Bigg) =0$$ let $u_1,u_2$ be two solutions. I want to show that $v := u_1 - u_2$ is a solution of the form $\mathrm{div} (A\nabla v)=0$ for a positive definite and symmetric matrix $A$.
The only advice given is to introduce the function $$ F(X):=\frac{X}{\sqrt{1+|X|^2}} $$ and use a green formula for $$ F(b)-F(a)=\int_0^1\partial_tF(ta+(1-t)b)\,\mathrm{d}t. $$
Attempt: When choosing $X=\nabla u_{1,2}$ and using linearity of the divergence I get
$$ \mathrm{div}\left (F(\nabla u_1)-F(\nabla u_2)\right ) = 0 = \mathrm{div}\int_0^1\partial_tF(t\nabla u_2+(1-t)\nabla u_1)\,\mathrm{d}t. $$
Here I am stuck as I don't find a suiting Green formula.
Now you need to compute. First note that $$ \partial_t F(t X + (1-t) Y) = DF(tX + (1-t)Y)(X-Y), $$ which in particular means that $$ \partial_t F(t \nabla u_2 + (1-t) \nabla u_1) = DF(t\nabla u_2 + (1-t)\nabla u_1)(\nabla u_2-\nabla u_1). $$
Returning to your formula, we note that $\nabla u_2 - \nabla u_1$ is independent of $t$, so
$$ 0 = \text{div} \left( \left[\int_0^1 DF(t\nabla u_2 + (1-t)\nabla u_1) dt\right] (\nabla u_2-\nabla u_1) \right) $$ and we get the desired result by setting $$ A = \int_0^1 DF(t\nabla u_2 + (1-t)\nabla u_1) dt. $$ To get more specific information about $A$ (ellipticity, etc), we need to compute $DF$. First note $$ \partial_j F(X)_k = \frac{ \sqrt{1+ |X|^2} \delta_{jk} - X_k X_j/ \sqrt{1 + |X|^2} }{1+ |X|^2}, $$ which we can rewrite as $$ DF(X) = \frac{I}{(1+|X|^2)^{1/2}} - \frac{X \otimes X}{(1+ |X|^2)^{3/2}}. $$ We can then plug in above to get more information about $A$.
EDIT: As requested, here is some more information about $A$.
It's obvious that $A$ is symmetric since $DF(\cdot)$ is. Let $z \in \mathbb{R}^n$ and note that $$ DF(X) z \cdot z = \frac{|z|^2}{(1+|X|^2)^{1/2}} - \frac{|X \cdot z|^2}{(1+ |X|^2)^{3/2}} = \frac{|z|^2 + |z|^2|X|^2 - |z\cdot X|^2 }{(1+ |X|^2)^{3/2}} \ge \frac{|z|^2}{(1+ |X|^2)^{3/2}}, $$ where the last inequality follows from Cauchy-Schwarz. Next note that $$ |t \nabla u_2 + (1-t) \nabla u_1| \le \Vert \nabla u_2\Vert_{L^\infty} + \Vert \nabla u_1\Vert_{L^\infty} =: C_0, $$ which means that $$ (1+ |t \nabla u_2 + (1-t) \nabla u_1|^2)^{3/2} \le (1+ C_0^2)^{3/2}. $$ Then $$ A z \cdot z \ge |z|^2 \int_0^1 \frac{dt}{(1+ |t \nabla u_2 + (1-t) \nabla u_1|^2)^{3/2}} \ge \frac{|z|^2}{(1+ C_0^2)^{3/2}}, $$ and we have that $A$ is elliptic with the ellipticity constant depending on $\nabla u_i$ via the constant $C_0$.