A nonlinear system with special structure 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?
 A: It seems unlikely that your problem is going to have a solution in the generality you've described but here goes.  There are two approaches you could try.
(1) Discretize the partial derivatives using finite differences on the grid.  At every point on the grid your PDE will give you a nonlinear equation.  You will have a total of $n^2$ equations in $n^2$ unknowns. But each equation will only contain a small number of unknowns. You can use a variant Newton's method to solve these equations.  But you have to work with the restrictions that you probably don't want to code up the Jacobian for the system, and you could not use a direct method to calculate the solution even if you did.  I would recommend looking at a Jacobian-Free Newton-Krylov method.
(2) Exploit the method of characteristics.  Rewrite your equation as 
$$
(cz+d) z_x - (az+b) z_y = 0.
$$
This gives characteristic equations 
$$
\frac{dx}{dt}=cz+d, \ \ \  \frac{dy}{dt}=a z+b, \ \ \ \frac{dz}{dt}= 0.
$$
There are numerical methods for solving nonlinear hyperbolic equations exploiting characteristics.  Sethian and Vladimirsky have a nice one.  Your problem does not quite fit into their scheme but their paper might help give you ideas.
So, if you problem does have a solution, one of these might work.  I would expect (1) to be more robust than (2), but also more expensive.
