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 (Boundary Condition): The ratio $\displaystyle\frac{z_x}{z_y}$ is known at the boundary. Thus, $z$ is known at the boundary.]

If we discretize $z_x$ and $z_y$, then we effectively have a large system of quadratic equations. Is there a good way to solve such systems?