Analytical solution for the following second order linear PDE I'm attempting to solve or to understand the conditions in which the following 1-dimensional second order, linear PDE can be solved. It is akin to a Fokker-Planck equation with constant coefficients, but with an additional linear sink term with space and time dependence:
$$ \frac{\partial f(x,t)}{\partial t} = - \mu \frac{\partial f(x,t)}{\partial x} + \frac{1}{2}\frac{\partial^2 f(x,t)}{\partial x^2} - u(x,t)f(x,t) $$
With natural boundary conditions in the domain $x\in(-\infty,\infty)$, $t\in[0,\infty)$, and initial condition 
$$ f(x,t=0) = \delta(x) $$
Additionally, $\mu$ is a positive real constant, and $u(x,t)$ is a strictly positive real function, monotonically increasing with $x$, and with asymptotes $a$ and $b$ at $x \rightarrow -\infty$, and $x \rightarrow \infty$ respectively (think of a sigmoid function that moves with time).
I have done an extense literature search that hasn't been fruitful. What I know is that the solution for the homogenous equation with $u(x,t)=0$ is
$$ f(x,t) = \frac{1}{\sqrt{2\pi t}} \exp\Bigl( -\frac{(x-\mu t)^2}{2t} \Bigr) $$
and that the solution with $u(x,t)=\alpha>0$ is
$$ f(x,t) = \frac{1}{\sqrt{2\pi t}} \exp\Bigl( -\frac{(x-\mu t)^2}{2t} - \alpha t\Bigr) $$
What is the best approach for solving this equation? Is there perhaps a change of variables that would simplify the problem? Is it even possible to obtain an analytical solution?

EDIT
I realized (also inspired by Harry49's answer) that the drift term can be removed with the following substitution:
$$ f(x,t) = g(x,t) \exp\Bigl( \mu x - \frac{\mu^2t}{2} \Bigr) $$
Which results in the following equation for $g$
$$ \frac{\partial g(x,t)}{\partial t} =  \frac{1}{2}\frac{\partial^2 g(x,t)}{\partial x^2} - u(x,t)g(x,t) $$
 A: Not a complete answer but I managed to convert it into an integral equation. I'm not sure if you can do better for a general $u$, but if you can, it will be essentially characterizing the semigroup generated by the operator $LG = -\frac{1}{2}\omega^2G-[U(\cdot,t)*G(\cdot,t)]$. Convolutions are pretty well-studied, but given that they are by definition global operators, I am sure the analysis is much harder and I am not familiar with the analysis needed to study such an operator in detail. I'm not sure how much can be done in closed form on that front, but I honestly doubt a nice analytic expression exists in general.
Starting from $$g_t = \frac{1}{2}g_{xx}-ug,$$ we take a Fourier transform $G:=\mathcal{F}g$ to obtain $$G_t = -\frac{1}{2}\omega^2G-[U(\cdot,t)*G(\cdot,t)](\omega)$$ with $G(\omega,0) = \mathcal{F}\delta = 1$. We can then use an integrating factor to obtain $$G(\omega,t) = e^{-\frac{1}{2}\omega^2t}-\int_{0}^te^{-\frac{1}{2}\omega^2(t-s)}[U(\cdot,s)*G(\cdot,s)](\omega)ds.$$
Now applying the inverse Fourier transform, we obtain $$g(x,t) = \frac{1}{\sqrt{2\pi t}}e^{-\frac{x^2}{2t}}-\int_0^t\mathcal{F}^{-1}\left\{e^{-\frac{1}{2}\star^2(t-s)}[U(\cdot,s)*G(\cdot,s)](\star)\right\}(x)ds.$$
We now focus on this inverse Fourier transform. We first split it into a product of the exponential and the convolution. Note that this convolution is in the Fourier variable. So 
$$
\mathcal{F}^{-1}\left\{e^{-\frac{1}{2}\star^2(t-s)}[U(\cdot,s)*G(\cdot,s)](\star)\right\}(x) = \left[\mathcal{F}^{-1}\left\{e^{-\frac{1}{2}\star^2(t-s)}\right\}(\diamond)*\mathcal{F}^{-1}\left\{[U(\cdot,s)*G(\cdot,s)](\star)\right\}(\diamond)\right](x).
$$
Notation is getting gross, but just keep in mind that the "inner" convolution was in the Fourier variable and the "outer" convolution is in the spatial variable. We can now compute these separate  inverse transforms.
$$
\mathcal{F}^{-1}\left\{e^{-\frac{1}{2}\star^2(t-s)}[U(\cdot,s)*G(\cdot,s)](\star)\right\}(x) = \left[\frac{1}{\sqrt{2\pi (t-s)}}e^{-\frac{\diamond^2}{2(t-s)}}*u(\diamond,s)g(\diamond,s)\right](x) = \int_{-\infty}^\infty\frac{1}{\sqrt{2\pi(t-s)}}e^{-\frac{(x-y)^2}{2(t-s)}}u(y,s)g(y,s)dy.
$$
Combining everything, we have
$$
g(x,t) = \frac{1}{\sqrt{2\pi t}}e^{-\frac{x^2}{2t}}-\int_0^t\frac{1}{\sqrt{2\pi(t-s)}}\left[\int_{-\infty}^\infty e^{-\frac{(x-y)^2}{2(t-s)}}u(y,s)g(y,s)dy\right]ds.
$$
This offloads a lot of the work to you, but to determine when the equation can be solved, you could examine under what conditions any of the steps I made in my derivation would not hold, like exchanging limits, assuming integrability, etc.
A: Let us introduce $g(X,t)$ such that
$$
g(X,t) = f(X+\mu t, t)\, e^{\int_0^t u(X+\mu\tau,\tau)\, \text d\tau} ,
$$
where $X=x-\mu t$.
Thus, we have the partial derivatives
$$
g_X = \left( f_X + f {\textstyle \int_0^t} u_X\, \text d\tau\right) e^{\int_0^t u\, \text d\tau} ,
$$
$$
g_{XX} = \left( f_{XX} + f {\textstyle \int_0^t} u_{XX}\, \text d\tau + 2 f_X {\textstyle \int_0^t} u_X\, \text d\tau + f\, ( {\textstyle \int_0^t} u_X\, \text d\tau)^2\right) e^{\int_0^t u\, \text d\tau} ,
$$
and 
$$
g_t = \left(\mu f_X + f_t + u f\right) e^{\int_0^t u\, \text d\tau} .
$$
Using the PDE, we have
$$
g_t -\tfrac12 g_{XX} = -\tfrac12 \left( f {\textstyle \int_0^t} u_{XX}\, \text d\tau + 2 f_X {\textstyle \int_0^t} u_X\, \text d\tau + f\, ( {\textstyle \int_0^t} u_X\, \text d\tau)^2\right) e^{\int_0^t u\, \text d\tau} .
$$
One recognizes the classical heat equation for $g$ if $u_X \equiv 0$, i.e. if $u$ is a function of $t$ only. Thus this method won't work in the general case where $u$ is a function of $(x,t)$.
For the general case, Laplace or Fourier transforms may provide integral representations of the solution (cf. [1] for the methodology).
It might be also of interests to look at the stationary solution.

[1] R. Habermann, Applied Partial Differential Equations with Fourier Series and Boundary Value Problems, 5th ed. Pearson Education Inc., 2013.
