I know how to use the characteristic method but I am in a little trouble with the following PDE: $$ \partial_xu(x,y)+\partial_yu(x,y)=a(x)u(x,y)u(x_0,y)\\ u(x,y=0)=f(x) $$ where $x_0$ is a fixed value belonging to the x domain and $f(x)$ is my boundary condition.
How could I apply the characteristic curves method in this case?
This is only a partial answer (as requested by OP in the comments to the question).
We start by treating $u(x_0,y)$ as a known function depending on $y$. Then the resulting linear first order PDE can be treated explicitely via method of characteristics resulting in: $$ u(x,y)=C(y-x)e^{\int_{x_0}^xa(s)u(x_0,s+y-x)ds} $$ for an arbitrary function $C$. Note that the lower bound for the intregation is arbitrary (different values are soaked up by the function $C$) but I went with $x_0$ because it will make the following steps easier.
Now basically you have two steps left. Two remove the implicit dependency of $u$ and to match the initial data $f$. As mentioned I got a partial result by tackling the implicit dependency of $u$ first:
If we evalute our solution at $x_0$ we achieve:
$$ u(x_0,y)=C(y-x_0)\cdot1 $$
Thus we get an explicit representation of $u(x_0,y)$ which we can plug back into the formula:
$$ u(x,y)=C(y-x)e^{\int_{x_0}^xa(s)C(s+y-x-x_0)ds} $$
If we evalute at $y=0$
$$ u(x,0)=C(-x)e^{\int_{x_0}^xa(s)C(s-x-x_0)ds}=f(x)\\ \Leftrightarrow \int_{z}^{-x_0}a(t+x0-z)C(t)dt=\ln(f(-z))-\ln(C(z)) $$
with $z=-x$ and $t=s-x-x_0$. The remaining task is to solve this integral equation for $C$, but I was only able to do this for constant $a$.
Assume $a$ is constant, then the equation reduces to $$ a\cdot\int_{z}^{-x_0}C(t)dt=\ln(f(-z))-\ln(C(z))\\ \Rightarrow -a\cdot C(z)=-\frac{f'(-z)}{f(-z)}-\frac{C'(z)}{C(z)} $$ This ODE can be solved explicitely and the solution reads: $$ C(z)=\frac{f(-z)}{-a\cdot\int_{-x_0}^z f(-s) ds+1} $$