Transformation not reproducible 
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.
 A: 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$.
