$$u_x + x \,y \,u_y + 2 x^2 \,z\,\ln z\,u_z = 0 $$
$$\frac{dx}{1}=\frac{dy}{xy}=\frac{dz}{2x^2z\ln|z|}=\frac{du}{0}$$
First family of characteristic curves from $\quad du=0 \quad\to\quad u=c_1$
Second family of characteristic curves from $\quad \frac{dx}{1}=\frac{dy}{xy} \quad\to\quad \frac{x^2}{2}-\ln|y|=c_2$
Third family of characteristic curves from $\quad\frac{dx}{1}=\frac{dz}{2x^2z\ln|z|} \quad\to\quad \frac{2}{3}x^3-\ln|\ln|z||=c_3 $
General solution of the PDE expressed on the form of implicit equation :
$$\Phi\left(u\:\:,\:\: \frac{x^2}{2}-\ln|y| \:\:,\:\: \frac{2}{3}x^3-\ln|\ln|z|| \right)=0$$
$\Phi$ is any differentiable function of three variables.
Or equivalently, general solution of the PDE on explicit form :
$$u(x,y,z)=\text{F}\left( \frac{x^2}{2}-\ln|y| \:\:,\:\: \frac{2}{3}x^3-\ln|\ln|z|| \right)$$
F is any differentiable function of two variables.
Although I have detailed a solution in the comments, I think it is probably good to provide the full solution here with an explanation.
Eliminating $dt$ in your system of ODEs yields the following form
$$\frac{dx}{x(y-u)} = \frac{dy}{y(u-x)} = \frac{du}{u(x-y)}$$
Noting that $dA/A = d \ln A$, this can be simplified to
$$\frac{d \ln x}{y-u} = \frac{d \ln y}{u-x} = \frac{d \ln u}{x-y}$$
Now, what the expression above is saying is that the vectors $(d \ln x, d \ln y, d \ln u)$ and $(y - u, u - x, x - y)$ are proportional to eachother i.e
$$(d \ln x, d \ln y, d \ln u) \propto (y - u, u - x, x - y)$$
You might notice that the RHS can be written as the curl of two vectors
$$(y - u, u - x, x - y) = (x, y, u) \times (1, 1, 1)$$
so what we really have is
$$(d \ln x, d \ln y, d \ln u) \propto (x, y, u) \times (1, 1, 1)$$
Remembering that the cross product of two vectors $a$ and $b$ yields a third vector $c$ which is orthogonal to both $a$ and $b$, the expression above really says that the tangent of the log of the solution curve is orthogonal to the vectors $(x, y, u)$ and $(1, 1, 1)$. So, if we take the dot product of the LHS with both vectors $(x, y, u)$ and $(1, 1, 1)$, we can get our integral curves
\begin{align}
(d \ln x, d \ln y, d \ln u) \cdot (x, y, u) &= ((x, y, u) \times (1, 1, 1)) \cdot (x, y, u) \\
&= 0 \\
\implies x d \ln x + y d \ln y + u d \ln u &= dx + dy + du \\
&= 0 \\
\implies x + y + u &= c_{1} \\\\
(d \ln x, d \ln y, d \ln u) \cdot (1, 1, 1) &= ((x, y, u) \times (1, 1, 1)) \cdot (1, 1, 1) \\
&= 0 \\
\implies d \ln x + d \ln y + d \ln u &= 0 \\
\implies \ln x y u &= c_{2} \\\\
\therefore \ln x y u &= f(c_{1}) \\
&= f(x + y + u)
\end{align}
and from here the result follows by applying the initial conditions as stated in the comments.
Note that this approach works well here because of the symmetry in the coefficients of the original PDE. Often times it can be quite difficult to determine the form of the two vectors whose curl is proportional to the tangent vector. However, I think this method gives a nice geometric view of what is happening.
Best Answer
Your calculus is correct. $$v\,v_x+v_y=-v\quad\text{is OK.}$$ The Charpit-Lagrange characteristic ODEs are : $$\frac{dx}{v}=\frac{dy}{1}=\frac{dv}{-v}$$ A first characteristic equation comes from solving $\frac{dy}{1}=\frac{dv}{-v}$ : $$v\,e^y=c_1 $$ A second characteristic equation comes from solving $\frac{dx}{v}=\frac{dv}{-v}$ : $$v+x=c_2$$ The general solution of the PDE on the form of implicit equation $c_1=f(c_2)$ is: $$v\,e^y=f(v+x)$$ $f$ is an arbitrary function. $$\boxed{v(x,y)=e^{-y}f\big(v(x,y)+x\big)}$$ $f$ has to be determined according to the boundary condition. $$u(x,0)=-x^2\quad\implies\quad v(x,0)=u_x(x,0)=-2x$$ $$-2x=e^0f(-2x+x)=f(-x)\quad\implies\quad f(X)=-2(-X)=2X$$ The function $f(X)$ is determined. We put it into the above general solution where $X=v+x$ : $$v=e^{-y}\:2(v+x)$$ $$v=\frac{2x}{e^y-2}$$ $$u=\int \frac{2x}{e^y-2} dx=\frac{x^2}{e^y-2}+C$$ $u(x,0)=-x^2\quad\implies\quad C=0$ $$\boxed{u(x,y)=\frac{x^2}{e^y-2}}$$