The general idea, as best I understand, is to take your differential identities, generate more by differentiating (prolongation), treat all derivatives as algebraically independent variables, and try to show that the putative identity is in the (algebraic) ideal generated by the known polynomials in these variables. The trick is to know how far to prolong.
While there are algorithmic approaches (I believe this goes by the name "Cartan Kuranishi algorithm, or something to that effect), I'll just assume that it will suffice to get all second derivatives of our defining equations. Here is the code.
The independent ariables and initial differential "polynomials".
ivars = {x, y};
dpolys = {D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] - Exp[-2*u[x, y]],
v[x, y] - (D[u[x, y], x, y] + D[u[x, y], x]*D[u[x, y], y])};
Get the Laplacian of v.
lapv = Total[Flatten[Outer[D[#1, {#2, 2}] &, {v[x, y]}, ivars]]];
Now get all first and second derivatives of these diff polys.
allpolys = Flatten[Join[dpolys, Outer[D, dpolys, ivars],
Outer[D[#1, {#2, 2}] &, dpolys, ivars], D[dpolys, x, y]]];
allvars = Variables[allpolys];
I want to impose an ordering on the "variables". I'm not sure this is necessary but at least it will help to make sure we use the same term ordering when we attempt to reduce the Laplacian of v to zero.
uvars = Select[allvars, ! FreeQ[#, u] &];
duvars = Select[uvars, ! FreeQ[#, Derivative] &];
vvars = Select[allvars, ! FreeQ[#, v] &]; dvvars =
Select[vvars, ! FreeQ[#, Derivative] &];
derivorder[Derivative[a : __][_][__]] := Plus[a]
sortuvars =
Append[Sort[duvars, derivorder[#1] >= derivorder[#2] &], u[x, y]];
sortvvars =
Append[Sort[dvvars, derivorder[#1] >= derivorder[#2] &], v[x, y]];
Small problem: when we go to compute a Groebner basis, the exponentials will not correctly be seen as invertible paramaters. So I replace them up front with a new symbol. We will treat it as invertible. It turns out that this both makes the reduction step viable and also makes the basis computation hugely faster.
allpolys2 = allpolys /. E^(-2 u[x, y]) -> exp2u;
Compute a(n algebraic) Groebner basis for our polynomials. Again, while they involve derivatives, we treat those as being algebraically independent variables.
Timing[
gb = GroebnerBasis[allpolys2, Join[sortvvars, sortuvars],
MonomialOrder -> DegreeReverseLexicographic,
CoefficientDomain -> RationalFunctions];]
(* Out[615]= {0.11, Null} *)
Now check whether the Laplacian of v lives in this (algebraic) ideal. If so, that proves that its vanishing is a consequence of the original differential equations.
PolynomialReduce[lapv, gb, Join[sortvvars, sortuvars],
MonomialOrder -> DegreeReverseLexicographic,
CoefficientDomain -> RationalFunctions][[2]]
(* Out[619]= 0 *)
[Last step: Declare victory and go home.]
v
and insert the explicit expression, but that gives a mass of symbols in which you'd need to replace $u_{xx}+u_{yy}$ by $\exp(-u)$ (not to mention the derivatives) and I don't see how to do that; I hope someone else will $\endgroup$