Suppose we have an $n \times n$ uniform grid, covering $[-1,1] \times [-1,1]$ (typically, n $\approx$ 500). We have a smooth, differentiable function $z(x,y)$ that we want to determine on the nodes of this grid. At every point, the function $z(x,y)$ satisfies a relation of the form:

$\displaystyle\frac{az + b}{cz + d} = \displaystyle\frac{z_x}{z_y}$

where $z_x$ and $z_y$ are the partial derivatives. Note that the values of $a, b, c, d$ are samples from four known smooth functions (so they are not constant at every point on the grid). How can we solve systems of equations with such structure?

**[EDIT 1]** (Boundary Condition): The ratio $\displaystyle\frac{z_x}{z_y}$ is known at the boundary. Thus, $z$ is known at the boundary.

**[EDIT 2]** Am I correct in understanding that this is a quasilinear first-order PDE and the method of characteristics will solve it up to level sets of z? Is there a robust way to solve such PDEs in the presence of noise and isolated singularities in a, b, c, d?