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 known, but 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]** Is this not a first-order, homogeneous, quasilinear PDE? Am I correct in assuming that the method of characteristics will solve this up to level sets of z?