Solve $y'(t) = a\left(e^{-k_1t}-e^{-k_2t}\right) - b\,y(t) - c\,y^2(t)$ I am attempting to resolve this first order nonlinear ODE
$$y'(t) = a\left(e^{-k_1t}-e^{-k_2t}\right) - b\,y(t) - c\,y^2(t)$$
where $a$, $b$, $c$, $k_1$ and $k_2$ are positive constants and $0\leq t<\infty$, with initial condition $y(0) = 0$.
Many moons ago, I typed the equation into Mathematica and I recall the solution being a ratio of differences of Bessel functions of different kinds. Vaguely from memory, it looked something like
$$y(t) = \frac{c_1K(\cdot)-c_2Y(\cdot)}{c_3J(\cdot)-c_4B(\cdot)}$$
where all capital letters are some form of Bessel functions, though I cannot remember the actual solution and I can't find where I wrote it down :(
For some reason, Mathematica no longer will solve this equation. So now I am attempting to solve the equation myself but have run into some roadblocks.
My Attempt
A priori the equation is a Riccati equation. So making the substitution
$$y=\frac{u'(t)}{c\,u(t)}$$
gives the following second order linear DE with variable coefficients
$$u''(t) + b\,u'(t) - ac\left(e^{-k_1t}-e^{-k_2t}\right)u(t)=0$$
My Mathematica does not like solving this one either. Power series isn't the best approach due to the last term in the equation. I can solve it using the Wronskian if I can find just one fundamental solution. Though I am not clever enough to find a fundamental solution to get the Wronskian method off the ground. However computation of the Wronskian is simple,
$$W_{u_1,u_2}(t)=Ce^{-bt}$$
where $C\neq0$ is arbitrary.
I then turned to Laplace transforms which looked promising at first
$$s^2U(s)-su(0)-u'(0)+b(sU(s)-u(0))-ac\left(U(s+k_1)-U(s+k_2)\right)=0$$
where $U(s)=\mathcal{L}_t\{u\}(s)$, but the shifts in the last two $U$ terms makes the solution mysterious to me.
Any suggestions on how to continue, or any alternative approaches to obtain the solution will be greatly appreciated.
 A: It may not be the best method, but one can at the very least write down a series solution to the original problem.
Let $$y(t) = \sum_{n=0}^{\infty} A_n t^n $$ 
Then $y'(t) = \sum_{n=0}^{\infty} (n+1)A_{n+1}t^n$, $y^2(t) = \sum_{n=0}^{\infty} \left( \sum_{k=0}^{n} A_k A_{n-k}\right) t^n$. So we arrive at
$$\sum_{n=0}^{\infty} (n+1)A_{n+1}t^n - a\sum_{n=0}^{\infty}\frac{(-1)^n(k_1^n - k_2^n)}{n!}t^n + b\sum_{n=0}^{\infty}A_nt^n+c\sum_{n=0}^{\infty} \left( \sum_{k=0}^{n} A_k A_{n-k}\right) t^n = 0$$ 
after substituting into the original problem. Gathering things up, and noting that $a_0 = 0$ from the initial condition, we have the following recurrence relation for the coefficients:
$$\begin{cases}
(n+1)A_{n+1} -a(-1)^n \frac{k_1^n-k_2^n}{n!} + bA_n + c \sum_{k=0}^{n} A_k A_{n-k} = 0 \\
A_0 = 0
\end{cases}$$
A couple of coefficients I found:
$A_1 = 0$
$A_2 = -\frac{1}{2}a(k_1 - k_2)$
$A_3 = \frac{1}{6}a(k_1 - k_2)(b + (k_1 + k_2))$
$A_4 = - \frac{1}{24}a(k_1 - k_2)(b(b+(k_1 + k_2))+k_1^2 + k_1 k_2 + k_2^2)$
EDIT: I wanted to make a quick mention for computing these coefficients efficiently since you're using Mathematica: use "memoization" to cut down on computation time, it allows you to compare the series solution to NDSolve pretty easily.
