My question is as follows: Let $A=\partial_x-\frac y2\partial_z$, $B=\partial_y+\frac x2\partial_z$, and $\Omega\subset \mathbb R^3$ be a smooth bounded open set. Take $g\in C^\infty(\Omega)$ (if you want even compactly supported in $\Omega$, although I don't think this is important).

Using the method of characteristics it is rather easy to solve the first order PDE $A(h)=g$, and in fact I obtained the solution (checked with mathematica) $$h(x,y,z)=\int_0^xg(\theta,y,\frac 12(xy+2z-\theta y))d\theta+H(y,\frac 12(xy+2z)).$$

Now, let $s\in C^\infty(\Omega)$ and assume $s$ to be very small in the $C^1(\Omega)$ norm, we can assume that $\|s\|_{C^1(\Omega)}\leq \varepsilon$, with $\varepsilon$ small enough. If it is of any use, we can assume that the $C^k(\Omega)$ norm of $\varepsilon$ is smaller than $\epsilon$, for $k\in\mathbb N$ big enough.

Now I would like to solve the perturbed PDE $$A(h)+B(s h)=g.$$ Again the method of characteristics works well if $s$ is a constant. With mathematica I found the explicit solution even if $s=s(x)$. For the general case however I doubt there is any hope of finding the explicit solution.

My question is: is there some general theory that ensures existence (I don't care about uniqueness) AND estimates of the solution $h$, in terms of $s$ and $g$? Even references are welcomed.

Many thanks

-Guido-