I have a PDE for function $f(x,y)$ in $(x,y) \in \mathbb{R}^2, x,y \geq 0$ where the PDE is $$\nabla f \cdot \vec{h} + f \nabla \cdot \vec{h} + \nabla^2 f = 0 \hspace{1em}$$ $$\vec{h} = \left(\begin{array}{c} -1 \\ 1 \\ \end{array}\right) \frac{1}{1 + a x - y}$$ with $a \in \mathbb{R}, 0 < a < 1$. The boundary conditions are imposed on $(x,y=0)$ and $(x=0, y)$. The expression $\nabla \cdot \vec{h} = \frac{1+a}{(1+a x - y)^2}$ is not zero.
I have not found this case in the book "Handbook of linear PDEs for engineers and scientists". I was thinking there may be a way to transform this PDE into a simpler Helmholtz equation. Do you know of any such transforms or other solution methods?
Update:
I have tried to find an approximate solution by considering $a x - y$ to be small. Thereby, the factor $\frac{1}{1 + a x - y}$ shoud be close to $1$. My approach was to replace the factor by $\frac{1}{1 + \eta(a x - y)}$ where I seek a solution in the form $$f = f_0(x,y) + \eta f_1(x,y) + \eta^2 f_2(x,y) + \dots $$ The solution of order $0$ is then determined by $$\left(a + \frac{1}{2} \right) g_0 + \nabla^2 g_0 = 0$$ when we substitute $f_{0} = g_{0} e^{\frac{1}{2} (x-y)}$.
However, is this approach admissible? After all, parameter $\eta$ has been introduced ad hoc. Is there a formal justification for this approach?