I trying to find the solution of the following Goursat problem for a second-order hyperbolic linear PDE $$ \begin{cases} u_{xy} + k(u_x+u_y) + (k^2 - \sigma^2 P(x-y))u = f(x,y) \\ u(x,0) = 0 \\ u(0,y) = 0 \end{cases} $$ where $P(x) = e^{-|x|/\tau}$, $k$, $\sigma$ and $\tau$ are constant, and $f(x,y)$ is a given function.
Do you have any ideas? Thanks for the help!
One way of tackling homogeneous equations of the form $u_{xy} + Au_x + Bu_y + Cu = 0$ is via the use of Laplace-Darboux transformations.
The technique is summarised in S.V. Meleshko Methods for Constructing Exact Solutions of Partial Differential Equations, Springer 2005.
A much more verbose exposition of the technique - with examples - can be found in Chapters XIII and XIV of A.R. Forsyth's Theory of Differential Equations, Part IV, Vol VI, Cambridge 1906.
After some fiddling I found the solution, which I'm posting for future reference.
First perform the transformation $$ u(x,y) = e^{-k(x+y)} v(x,y) $$ which essentially sets $k=0$ by giving the simpler Goursat problem for $v(x,y)$ $$ \begin{cases} v_{xy} - \sigma^2 P(x-y) v = e^{k(x+y)}f(x,y) \\ v(x,0) = 0 \\ v(0,y) = 0 \end{cases} \tag{1}\label{1} $$
At this point one can use the Riemann method to find the solution to the inhomogeneous problem. In this case it amounts to finding the solution $A(x_0,y_0;x,y)$ to the homogeneous PDE $$ \begin{cases} A_{x_0y_0} - \sigma^2 P(x_0-y_0)A = 0 \\ A(x_0,y;x,y) = 1 \\ A(x,y_0;x,y) = 1 \end{cases} \tag{2}\label{2} $$
Suppose at this point that there was no absolute value in $P(x)$. In this case the PDE would be $$ \begin{cases} A_{x_0y_0} - \sigma^2 e^{-x_0/\tau} e^{y_0/\tau} A = 0 \\ A(x_0,y;x,y) = 1 \\ A(x,y_0;x,y) = 1 \end{cases} \tag{3}\label{3} $$ By performing the change of variables $z=e^{-x_0/\tau}$, $w=e^{y_0/\tau}$ and using the known results of the PDE $A_{zw} + \lambda A =0$ one can check that the solution of the Goursat problem \eqref{3} is $$ A_0(x_0,y_0;x,y) = J_0\left(2\sigma\tau\sqrt{\left(e^{-x_0/\tau}-e^{-x/\tau}\right)\left(e^{y_0/\tau}-e^{y/\tau}\right)}\right) $$ The solution if the absolute value is replaced with the opposite of the argument is by similar reasoning $A_0(-x_0,-y_0,-x,-y)$. It is clear then that the solution of \eqref{2} is $$ A(x_0,y_0;x,y) = \theta(x_0-y_0)A_0(x_0,y_0;x,y) + \theta(y_0-x_0) A_0(-x_0,-y_0,-x,-y) $$
By the Riemann's formula the solution to \eqref{1} is $$ v(x,y) = \int_0^x dx_0 \int_0^y dy_0 e^{k(x_0+y_0)}f(x_0,y_0)A(x_0,y_0;x,y) $$ so the final answer to the original question is $$ u(x,y) = e^{-k(x+y)}\int_0^x dx_0 \int_0^y dy_0 e^{k(x_0+y_0)}f(x_0,y_0)A(x_0,y_0;x,y) $$
