How do I solve a PDE with multiple Dirac functions? I am exposed to a PDE in the following form: 
$\frac{\partial f}{\partial t}=\alpha \frac{\partial^2 f}{\partial x^2}-\beta \frac{\partial f}{\partial x} + \mu_a P_a(t) \delta(x-1)+ \mu_b P_b(t) \delta(x-N+1)$
I'm trying to find two linearly independent solutions. I know that the time derivative could be treated with Laplace transform and that I should divide the solution interval according to dirac arguments. Although, I don't get the expected result.  
Could someone help me to proceed in solving this?
It is about applying Feller results on a modified Fokker-planck (Chapman-Kolmogovo) which has the same form of the given equation.  
Edit: 
$P_a$ and $P_b$ are defined by their differential equations as follows: 
 \begin{equation}
 \frac{dP_a}{dt}=-\mu_a P_a + \lim _{x\rightarrow 0} \left( \alpha \frac{\partial^2 f}{\partial x^2}-\beta \frac{\partial f}{\partial x}\right)
 \end{equation}
 \begin{equation}
\frac{dP_b}{dt}=-\mu_b P_b + \lim _{x\rightarrow N} \left( \alpha \frac{\partial^2 f}{\partial x^2}-\beta \frac{\partial f}{\partial x}\right)
 \end{equation}
Actually I am only interested in fundamental solutions of the PDE, I am not sure if the presence of dirac function affects the two fundamental solutions $e^{\xi_i}$ given that $\xi_i$ are roots of the associated characteristic equation. 
 A: $\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}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\Li}[1]{\,\mathrm{Li}_{#1}}
 \newcommand{\ol}[1]{\overline{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

$\ds{\partiald{\,\mathrm{f}\pars{x,t}}{t} =
\alpha\,\partiald[2]{\,\mathrm{f}\pars{x,t}}{x} - \beta\,\partiald{\,\mathrm{f}\pars{x,t}}{x}\ +\mu_{a}\,\mathrm{P}_{a}\pars{t}\delta\pars{x - 1} + \mu_{b}\,\mathrm{P}_{b}\pars{t}\delta\pars{x - N + 1}}$.

Note that 
\begin{align}
&\pars{\partiald{}{t} - {\beta^{2} \over 4\alpha}}
\,\mathrm{f}\pars{x,t}
=
\alpha\,\pars{%
\partiald{}{x} - {\beta \over 2\alpha}}^{2}\,\mathrm{f}\pars{x,t} +
\mu_{a}\,\mathrm{P}_{a}\pars{t}\delta\pars{x - 1} + \mu_{b}\,\mathrm{P}_{b}\pars{t}\delta\pars{x - N + 1}\end{align}

With $\ds{\,\mathrm{F}\pars{x,t} \equiv \exp\pars{-\,{\beta \over 2\alpha}\,x - {\beta^{2} \over 4\alpha}\,t}\,\mathrm{f}\pars{x,t}}$:
\begin{align}
\partiald{\,\mathrm{F}\pars{x,t}}{t} & = \alpha\,\partiald[2]{\,\mathrm{F}\pars{x,t}}{x} + \,\mathrm{g}_{a}\pars{t}\delta\pars{x - 1} + \,\mathrm{g}_{b}\pars{t}\delta\pars{x - N + 1}\tag{1}
\\[4mm] & \mbox{where}\
\left\lbrace\begin{array}{rcl}
\ds{\,\mathrm{g}_{a}\pars{t}} & \ds{=} &
\ds{\mu_{a}\exp\pars{-\,{\beta \over 2\alpha}}\exp\pars{-\,{\beta^{2} \over 4\alpha}\,t}}\,\mathrm{P}_{a}\pars{t}
\\[2mm]
\ds{\,\mathrm{g}_{a}\pars{t}} & \ds{=} &
\ds{\mu_{b}\exp\pars{-\,{\beta \over 2\alpha}\,\bracks{N - 1}}\exp\pars{-\,{\beta^{2} \over 4\alpha}\,t}}\,\mathrm{P}_{b}\pars{t}
\end{array}\right.
\end{align}

Now, I'll 'take' Laplace Transform in both sides of $\pars{1}$:
\begin{align}
-\,\mathrm{F}\pars{x,0} + s\,\hat{\,\mathrm{F}}\pars{x,s} & =
\alpha\,\partiald[2]{\hat{\,\mathrm{F}}\pars{x,s}}{x} +
\hat{\,\mathrm{g}}_{a}\pars{s}\delta\pars{x - 1} + \hat{\,\mathrm{g}}_{b}\pars{s}\delta\pars{x - N + 1}
\end{align}
which leads to
\begin{align}
\pars{\partiald[2]{}{x} - {s \over \alpha}}\hat{\,\mathrm{F}}\pars{x,s}
 & =
-\,{1 \over \alpha}\bracks{\,\mathrm{F}\pars{x,0} +
\hat{\,\mathrm{g}}_{a}\pars{s}\delta\pars{x - 1} + \hat{\,\mathrm{g}}_{b}\pars{s}\delta\pars{x - N + 1}}
\end{align}

In terms of the Green's Function $\ds{\,\mathrm{G}\pars{s,x,x'}}$ the solution,
for $\ds{\hat{\,\mathrm{F}}\pars{x,s}}$, is written as
\begin{align}
\hat{\,\mathrm{F}}\pars{x,s} & =
\varphi\pars{x,s}
\\[3mm] & -\,{1 \over \alpha}\int_{-\infty}^{\infty}\,\mathrm{G}\pars{s,x,x'}\bracks{\,\mathrm{F}\pars{x',0} +
\hat{\,\mathrm{g}}_{a}\pars{s}\delta\pars{x' - 1} + \hat{\,\mathrm{g}}_{b}\pars{s}\delta\pars{x' - N + 1}}\,\dd x'
\\[8mm] & =
\varphi\pars{x,s} -
{1 \over \alpha}\bracks{\hat{\,\mathrm{g}}_{a}\pars{s}\,\mathrm{G}\pars{s,x,1} + \hat{\,\mathrm{g}}_{b}\pars{s}\,\mathrm{G}\pars{s,x,N - 1}}
\\[4mm] & -
{1 \over \alpha}\int_{-\infty}^{\infty}\,\mathrm{G}\pars{s,x,x'}
\,\mathrm{F}\pars{x',0}\,\dd x'
\end{align}
$\ds{\varphi\pars{x,s}}$ satisfies
$\ds{\pars{\partiald[2]{}{x} - {s \over \alpha}}\varphi\pars{x,s} = 0}$ and the
$\ds{x}$-boundary conditions of $\ds{\hat{\,\mathrm{F}}\pars{x,s}}$ ( lets
assume, for example, that it occurs at $\ds{x = 0}$ ).
$$
\pars{\partiald[2]{}{x} - {s \over \alpha}}{\,\mathrm{G}}\pars{s,x,x'} = \delta\pars{x - x'}\,,\qquad\,\mathrm{G}\pars{s,0,x'} = 0
$$
The general solution is given by
$$
\,\mathrm{G}\pars{s,x,x'} =
\left\lbrace\begin{array}{lcl}
\ds{0} & \mbox{if} & \ds{x < x'}
\\[2mm]
\ds{-\root{\alpha \over s}\sinh\pars{\root{s \over \alpha}\bracks{x - x'}}}
& \mbox{if} & \ds{x > x'}
\end{array}\right.
$$
A: I am not sure whether you will find this helpful, but I would attack this problem as follows.
First step
First, I would try to get rid of the differential operators in $x$ by performing a partial Fourier transform, i.e. a Fourier transform only in $x$ and not in $t$. This would convert these differential operators into algebraic monomials.
Concretely, if we use the tilde to denote the partial Fourier transform
$$\tilde f (t,y) = \frac 1 {\sqrt {2 \pi}} \int \limits _{-\infty} ^\infty \Bbb e ^{- \Bbb i y x} f(t,x) \ \Bbb d x$$
(in the sense of distributions), then
$$\widetilde {\frac {\partial f} {\partial x}} (t,y) = - \Bbb i y \tilde f (t,y), \quad \widetilde {\frac {\partial^2 f} {\partial x^2}} (t,y) = - y^2 \tilde f (t,y), \quad \widetilde {\delta_a} (y) = \frac 1 {\sqrt {2 \pi}} \Bbb e ^{- \Bbb i y a}$$
(where $\delta_a (x) = \delta (x-a)$), so (using that $\widetilde {\frac {\partial f} {\partial t}} = \frac {\partial \tilde f} {\partial t}$) this Fourier transform would transform the initial equation into
$$\frac {\partial \tilde f} {\partial t} (t,y) = (- \alpha y^2 + \beta \Bbb i y) \tilde f (t,y) + \frac 1 {\sqrt {2 \pi}} \mu_a P_a (t) \Bbb e ^{- \Bbb i y} + \frac 1 {\sqrt {2 \pi}} \mu_b P_b (t) \Bbb e ^{- \Bbb i y (N-1)}$$
which is an inhomogeneous differential equation of order $1$ that you can solve with the standard tool.
Second step
To solve the above equation, first treat the homogeneous case, thus considering the equation
$$\frac {\partial \tilde f} {\partial t} (t,y) = (- \alpha y^2 + \beta \Bbb i y) \tilde f (t,y) .$$
This has the solution
$$\tilde f (t,y) = C \Bbb e ^{(- \alpha y^2 + \beta \Bbb i y) t}$$
where $C>0$ is some arbitrary constant.
Next, treat the inhomogeneous equation by considering $C$ to be a function $C(t)$, and not a constant anymore, so
$$\tilde f (t,y) = C(t) \Bbb e ^{(- \alpha y^2 + \beta \Bbb i y) t} .$$
Replacing this formula for $\tilde f$ into the inhomogeneous equation leads to the new differential equation
$$C'(t) = \Bbb e ^{- (- \alpha y^2 + \beta \Bbb i y) t} \left( \frac 1 {\sqrt {2 \pi}} \mu_a P_a (t) \Bbb e ^{- \Bbb i y} + \frac 1 {\sqrt {2 \pi}} \mu_b P_b (t) \Bbb e ^{- \Bbb i y (N-1)} \right)$$
which can be integrated directly to yield
$$C(t) = \int \limits _{t_0} ^t \Bbb e ^{- (- \alpha y^2 + \beta \Bbb i y) s} \left( \frac 1 {\sqrt {2 \pi}} \mu_a P_a (s) \Bbb e ^{- \Bbb i y} + \frac 1 {\sqrt {2 \pi}} \mu_b P_b (s) \Bbb e ^{- \Bbb i y (N-1)} \right) \Bbb d s$$
with $t_0$ arbitrary in the correct domain of definition required by $f$ (that you do not give, so I shall assume $f$ to be defined on $\Bbb R^2$ - in which case $t_0 \in \Bbb R$). Plugging this back in the formula of $\tilde f (t,y)$ yields
$$\tilde f (t,y) = \Bbb e ^{(- \alpha y^2 + \beta \Bbb i y) t} \ \int \limits _{t_0} ^t \Bbb e ^{- (- \alpha y^2 + \beta \Bbb i y) s} \left( \frac 1 {\sqrt {2 \pi}} \mu_a P_a (s) \Bbb e ^{- \Bbb i y} + \frac 1 {\sqrt {2 \pi}} \mu_b P_b (s) \Bbb e ^{- \Bbb i y (N-1)} \right) \Bbb d s .$$
Third step
To obtain $f$, perform an inverse Fourier transform in the $y$ variable, obtaining
$$f(t,x) = \frac 1 {\sqrt {2 \pi}} \int \limits _{-\infty} ^\infty \Bbb e ^{\Bbb i x y} \tilde f (t,y) \ \Bbb d y = \\
\frac 1 {2 \pi} \int \limits _{-\infty} ^\infty \int \limits _{t_0} ^t e ^{\Bbb i x y} \Bbb e ^{(- \alpha y^2 + \beta \Bbb i y) t} \ \Bbb e ^{- (- \alpha y^2 + \beta \Bbb i y) s} \left( \mu_a P_a (s) \Bbb e ^{- \Bbb i y} + \mu_b P_b (s) \Bbb e ^{- \Bbb i y (N-1)} \right) \Bbb d s \ \Bbb d y.$$
Maybe this can be further slightly simplified; I do not know if it answers your question in the way you were hoping for, but if not it can at least suggest an approach. Not having concrete formulae for $P_a$ and $P_b$ necessarily leaves things in a somewhat unfinished, generic state.
