I need to solve the following ODE (a steady state Fokker-Planck/Kolmogorov Forward equation)
\begin{equation} \dfrac{\partial^2}{\partial x^2}\left[ax^2f(x)\right] - \dfrac{\partial}{\partial x}\left[bxf(x)\right] - c f(x) + c\delta(x-\hat{x}) = 0 \end{equation}
where $x\in(0,\infty)$ and where $f(x)\rightarrow 0$ as $x \rightarrow \infty$. The constants $a,b,c$ are all positive and $\hat{x}$ is positive and far away from zero (far enough such that the solution at $x=0$ should be $f(0)=0$). The solution ignoring the Dirac impulse $\delta(x-\hat{x})$ (the solutions where $x \neq \hat{x}$) and imposing that $\lim_{x \rightarrow \infty}f(x) = \lim_{x \rightarrow 0}f(x) = 0$ is given by
$$ f(x) = \begin{cases} c_1 x^{\lambda_1} \quad \text{if} \ x<\hat{x}\\ c_2 x^{\lambda_2} \quad \text{if} \ x>\hat{x}. \end{cases} $$
This is found by guessing the solution has the shape $Ax^\lambda$ and plugging such guess to the ODE shown above. Doing this we can verify that the solution satisfies the ODE and gives us the expression for $\lambda$.
$$ \lambda = \dfrac{b-3a \pm \sqrt{(b-a)^2+4ac}}{2a} $$
Since I have a Dirac impulse at $x=\hat{x}$, I should be solving for two ODEs, one below $x=\hat{x}$ and another above $x=\hat{x}$. I then have to glue together both solutions at $x=\hat{x}$. This tells us that
$$c_1 = c_2 \hat{x}^{\lambda_2-\lambda_1}.$$
Finally, as we have a Dirac delta at $x=\hat{x}$, we get yet another boundary condition by integrating the ODE around $\hat{x}$.
$$ \lim_{\epsilon \rightarrow 0}\int_{\hat{x}-\epsilon}^{\hat{x}+\epsilon} \left(\dfrac{\partial^2}{\partial x^2}\left[ax^2f(x)\right] - \dfrac{\partial}{\partial x}\left[bxf(x)\right] - c f(x) + c\delta(x-\hat{x}) \right) \text{d}x = 0$$
When we do this, the second and third terms of the ODE drop out (by continuity these two terms must be equal to zero when we integrate them around $\hat{x}$). Thus we end up with the next condition.
$\dfrac{c}{a} + c_2\lambda_2\hat{x}^{\lambda_2+1} = c_1\lambda_1\hat{x}^{\lambda_1+1}$
Putting everything together we get the coefficients $c_1,c_2$
$$ c_1 = \dfrac{c}{a}\dfrac{\hat{x}^{-(1+\lambda_1)}}{(\lambda_1-\lambda_2)}, \ c_2 = \dfrac{c}{a}\dfrac{\hat{x}^{-(1+\lambda_2)}}{(\lambda_1-\lambda_2)},$$
and we use the $\lambda$ that is a negative root for $x > \hat{x}$, (and viceversa for $x<\hat{x}$).
The steps I followed to tackle this question are based on two previous questions I asked here. One case did not have diffusion but had convection (just a first order ODE - link here), whereas the other case had diffusion but no convection (second order ODE without the first derivative term - link here). This case now has diffusion and convection.
The problem I have is that 1) The solution does not match comparisons with numerical results which I know are correct. 2) The solution does not integrate to 1, that is $\int f(x)\text{d}x \neq 1$. I tried to impose this as a boundary condition but once again my results do not make much sense. Where did I make a mistake, or how would you solve this problem?
A bit more on the intuition behind the problem. The ODE in question is a steady state Fokker-Planck (or Kolmogorov Forward) Equation. Mass is injected at $x=\hat{x}$ and dissipates both to the left and right of $x=\hat{x}$ because of diffusion. Diffusion is controled by $a$. The result should display a density with mass skewed to the right of $\hat{x}$ since $b>0$. $b$ controls convection. Then, mass anywhere in $x\in(0,\infty)$ is taken out at a rate $c$ and reinjected back to $x=\hat{x}$.