Solution to two-dimensional PDE (wave/Klein-Gordon type equation)

Question

I was playing with an optimization problem and ended up reducing it to solving the following PDE: $$ a^2 xy \frac{\partial^2 f}{\partial x \partial y} + f = -axy(y+2ax) $$ where $a \in \mathbb{R}^+$ is fixed and $f : \mathbb{R}^2 \to \mathbb{R}$ is unknown. I have no experience whatsoever with PDEs, but noticed I could separate $x$ and $y$ by changing variables $t = (x+y)/2$ and $r = (x-y)/2$, giving $$ \left(\frac{\partial^2}{\partial^2 t} - \frac{\partial^2}{\partial^2 r} + \frac{4}{a^2(t^2-r^2)} \right) f = g(t,r) $$ for some fixed $g$ (not relevant here). This looks something like the 1+1-dimensional (time + 1D) Klein-Gordon equation $$ \left(\frac{\partial^2}{\partial^2 t} - \frac{\partial^2}{\partial^2 r} + m^2 \right) f = 0 $$ except for the rightmost term in the bracket and the RHS of the equation, which both vary with $t,r$ in my case. It also looks a little bit like the inohomogeneous Helmholtz equation or the screened Poisson equation depending on whether $xy > 0$ or $xy < 0$, except they also have a constant rightmost term in the bracket AND have no time-component. I have no idea where to start, having no expertise in PDEs. I tried looking at the following papers but had a hard time following. Any guidance, hints, solutions most appreciated.