Find the common solution of $$ u_{xy}+a(x,y)u_x=0~~~~~\text{in}~~~~~\Omega:=\left\{(x,y)\in\mathbb{R}^2 : \lvert x-x_0\rvert<d_1, \lvert y-y_0\rvert <d_2\right\},\\ d_1,d_2>0, (x_0,y_0)\in\mathbb{R}^2, a\in C(\Omega) $$
In the meantime, I thought about an answer to this question (see below). Would be great to get a feedback to it.
With kind regards
I add my answer in order to get a feedback from you. Please tell me, if my general solution is right or maybe nonsense, would be very kind of you!
First of all I use the theorem of Schwarz and substitute $z:=u_x$, getting an ODE of order 1 in $y$ (here $x$ is only a parameter): $$ z_y=-a(x,y)z~~~~~(1) $$ If one considers $a$ as only dependent on $y$ and $x$ only as a parameter, then $a$ is continious on the intervall $I:=(y_0-d_2,y_0+d_2)$. So the Fundamental Theorem of Calculus says that $$ A(x):=\int_{y_0}^{y}a(x,\tau)\ d\tau~~~~~(2) $$ is an antiderivative on $I$. So I continue (1) by separation of variables, getting $$ z=C_1(x)\cdot\exp\left(-\int_{y_0}^{y}a(x,\tau)\ d\tau\right). $$ Resubstituting, it follows $$ u(x,y)=\int_{x_0}^{x}C_1(\eta)\exp\left(-\int_{y_0}^y a(\eta,\tau)\, d\tau\right)\, d\eta+C_2(y)~~~~~(3) $$ with $C_1\in C^1((x_0-d_1,x_0+d_1))$ and $C_2\in C^1(I)$ arbitrary functions.
Comment: Instead of chosing $y_0$ and $x_0$ as lower borders in (2) and (3)it is of course possible to chose any point in $I$ respectively any point in $(x_0-d_1,x_0+d_1)$.
$\newcommand{\angles}[1]{\left\langle #1 \right\rangle}% \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}% \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}% \newcommand{\dd}{{\rm d}}% \newcommand{\isdiv}{\,\left.\right\vert\,}% \newcommand{\ds}[1]{\displaystyle{#1}}% \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}% \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}% \newcommand{\ic}{{\rm i}}% \newcommand{\imp}{\Longrightarrow}% \newcommand{\pars}[1]{\left( #1 \right)}% \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}}% \newcommand{\sgn}{\,{\rm sgn}}% \newcommand{\ul}[1]{\underline{#1}}% \newcommand{\verts}[1]{\left\vert #1 \right\vert}% \newcommand{\yy}{\Longleftrightarrow}$
\begin{align} \partiald{{\rm u}\pars{x,y}}{x} &= \exp\pars{-\int_{y_{0}}^{y}{\rm a}\pars{x,y'}\,\dd y'} + \partiald{{\rm u}\pars{x,y_{0}}}{x} \\ \partiald{{\rm u}\pars{x,y_{0}}}{x} &\equiv \left.\partiald{{\rm u}\pars{x,y'}}{x}\right\vert_{y' = y_{0}} \equiv \tilde{\rm u}\pars{x,y_{0}} \end{align}
\begin{equation} {\rm u}\pars{x,y} = \int_{x_{0}}^{x} \exp\pars{-\int_{y_{0}}^{y}{\rm a}\pars{x',y'}\,\dd y'}\,\dd x' + \int_{x_{0}}^{x}\tilde{\rm u}\pars{x',y_{0}}\,\dd x' + {\rm u}\pars{x_{0},y} \tag{1} \end{equation}
In addition, $$ {\rm u}_{x}\pars{x,y}{\rm u}_{xy}\pars{x,y} + {\rm a}\pars{x,y}{\rm u}_{x}^{2}\pars{x,y}=0\,, \quad\imp\quad \bracks{{1 \over 2}\,\partiald{}{y} + {\rm a}\pars{x,y}}{\rm u}_{x}^{2}\pars{x,y} = 0 $$$$ \partiald{}{y}\bracks{% \exp\pars{2\int_{y_{0}}^{y}{\rm a}\pars{x,y'},\dd y'} {\rm u}{x}^{2}\pars{x,y}} = 0 $$ and $$ {\rm u}{x}^{2}\pars{x,y} = \tilde{\rm u}^{2}\pars{x,y_{0}}\exp\pars{-2\int_{y_{0}}^{y}{\rm a}\pars{x,y'},\dd y'} $$
Let's assume we have two solutions $\phi\pars{x,y}$ and $\varphi\pars{x,y}$ of the differential equation for given values of ${\rm u}\pars{x_{0},y}$ and $\tilde{\rm u}\pars{x,y_{0}}$. Then, we have $\delta_{x}^{2}\pars{x,y} = 0$ where $\delta\pars{x,y} \equiv \phi\pars{x,y} - \varphi\pars{x,y}$. We conclude that $$ \delta\pars{x,y} = {\rm f}\pars{y}\ \mbox{where}\ {\rm f}\ \mbox{is an }\ {\it arbitrary}\ \mbox{function}. $$ However, $$ {\rm f}\pars{y} = \phi\pars{x_{0},y} - \varphi\pars{x_{0},y} = {\rm u}\pars{x_{0},y} - {\rm u}\pars{x_{0},y} = 0 \quad\imp\quad \phi\pars{x,y} = \varphi\pars{x,y} $$
We conclude that the solution $\pars{1}$ is the unique solution for a given boundary conditions as defined by $\ds{{\rm u}\pars{x_{0},y}}$ and $\ds{\tilde{\rm u}\pars{x_{0},y} \equiv \left.\partiald{{\rm u}\pars{x,y'}}{y'}\right\vert_{y'\ =\ y_{0}}}$.
