Si sólo estás interesado en cómo resolver este tipo de ecuaciones, puede que encuentres esto útil.
No hay necesidad de utilizar este tipo de cambio de variables "de sopetón" si te introduces en el método de las características para resolver EDP de primer orden. En tu caso, la EDP también es lineal y puede escribirse como sigue:
$$a(x,y) \, p + b(x,y)\, q = f(x,y,z), \quad p = \partial_x z, \quad q = \partial_y z,$$ y $a = 1, b = 1$ son funciones conocidas de sus variables independientes y $f = 0$ es una función lineal sobre $z$ . El método de las características dice entonces:
$$\frac{\mathrm{dx}}{a} = \frac{\mathrm{dy}}{b} = \frac{\mathrm{dz}}{f} \Rightarrow \frac{\mathrm{dx}}{1} = \frac{\mathrm{dy}}{1} = \frac{\mathrm{d}z}{0}.$$ Las dos primeras igualdades nos lo dicen:
$$ \mathrm{d}x - \mathrm{d}y = 0,$$ así que $x-y=c$ (constante) en algún lugar . La fracción $\mathrm{d}z/0$ lleva a establecer $z=k$ (constante), de nuevo, en algún lugar . Si pones $k$ en función de $c$ Entonces la solución se convierte en:
$$z = k(c) = k(x-y),$$
donde $k(\cdot)$ es una función arbitraria de $x-y$ es decir, la curva característica de la ecuación.
Espero que esto ayude.
Salud.