I'm interested in finding a solution (or an approximate solution) to the following PDE:
$f_t + (a-bt^2)f_v = c\delta(v)$
This is the Vlasov equation for a distribution function $f$ with no spatial dependence, a delta function source term, and an acceleration of the form $(a-bt^2)$.
Formally speaking, one can always solve an equation like this. Just apply the method of characteristics to the general equation
$$f_t+(a-bt^2)f_v=g(t,v)$$ which reads $$dt=\frac{dv}{a-bt^2}=\frac{df}{g(t,v)}$$ The next step is to find two distinct integrals of motion. These are readily found since
$$dv-(a-bt^2)dt=0\iff v-at+\frac{bt^3}{3}=C_1$$
and with the aid of this integral of motion we can find a second one:
$$df=g(t,v)dt\iff f-\int^tdt'g(t',v-at+\frac{bt^3}{3}+at'-\frac{bt'^3}{3})=C_2$$
From these then one can find the most general solution to the equation above to be
$$f(t,v)=\int^tdt'g(t',v-at+\frac{bt^3}{3}+at'-\frac{bt'^3}{3})+ H\left(v-at+\frac{bt^3}{3}\right)$$
where $H$ is an arbitrary function. The first term is a solution to the inhomogeneous equation and the 2nd one is a solution to it's homogeneous part. To solve your particular problem, $g(t,v)=c\delta(v)$ you need to solve a cubic equation, which is not fun because of the casework and the lack of physical interpretation and boundary conditions. Regardless, the form of the solutions is going to look like
$$f(t,v)=\frac{c}{a}\sum\frac{\theta\left(t-2\sqrt{\frac{a}{b}}r_i(X)\right)}{|1-4r_i^2(X)|}$$ where $r_i(x)$ represents the set of real roots of the equation $4u^3-3u-x=0$ and $X$ being the following dimensionless quantity
$$X=\frac{3b^{1/2}[v-at+b\frac{t^3}{3}]}{2a^{3/2}}$$
EDIT: Now that the physical picture has been added it is much easier to provide a concrete solution based on the above. We posit that at time zero, there are no particles present, thus $f(0,v)=0$. Applying the boundary condition to the general solution we find after a short calculation that the solution is
$$f(t,v)=c\int_0^tdt'\delta(v-at+\frac{bt^3}{3}+at'-\frac{bt'^3}{3})$$
which yields a solution exactly as above but with the added restriction that the roots of the equation have to be taken into account are the ones that are positive!
This picture is supported further when one solves the problem by physical considerations only, monitoring the velocity of every particle generated at an arbitrary time by the source. The velocity of the particle generated at time $T$ as a function of time is given by
$$v(t,T)=a(t-T)-\frac{b}{3}(t^3-T^3)$$
We ask: At a certain time t, given that the stationary particle production rate is $\frac{dN}{dT}=c$, how many particles are contained in the velocity interval $(v,v+dv)$? It obviously is the ones that were generated at the correct moments $(T, T+dT)$ so that their velocities are the ones required. Formally
$$dN(v\leq v(t,T)\leq v+dv)=dN(v-at+\frac{bt^3}{3}\leq -aT+\frac{bT^3}{3}\leq v-at+\frac{bt^3}{3}+dv)$$
which yields that $T$ must lie in the (possibly multiple) intervals $r\leq T\leq r+\frac{dv}{br^2-a}$ where $r$ is a positive (since time of production is positive) root of the equation $-ar+\frac{br^3}{3}=v-at+\frac{bt^3}{3}$. Importantly, this root also has to be less than $t$ itself, since the time of production of the particle has to precede the time at which a particle is measured to have the purported velocity $v$.With these in mind, from the rate of particle production we get that $$dN(v\leq v(t,T)\leq v+dv)=\sum_n dN(r_n\leq T\leq r_n+\frac{dv}{br_n^2-a})=cdv\sum_n\frac{1}{|br_n^2-a|} \\~,~r_n\in[0,t]$$
which finally implies
$$\frac{dN}{dv}=f(t,v)=c\sum_n\frac{1}{|br_n^2-a|}$$
This is precisely what we expect from the integration of the delta functions of each distinct positive root of the cubic presented before the edit (modulo rescaling of various prefactors).
Of course, to obtain the global properties of this distribution we need to do more work. In particular, we need to monitor how much time it actually takes for the first particle to arrive at a particular velocity. Before that time the distribution is actually identically zero at that point. Assuming for simplicity that $a,b>0$ and sparing you the details we have that for $v<0 ~, ~t\geq t_0(v)$ where $t_0$ is the only positive solution of the equation $at_0-\frac{bt_0^3}{3}=v$, for $v>\frac{2a}{3}\sqrt{\frac{a}{b}}$, the distribution is identically zero for all times, and for $0\leq v\leq\frac{2a}{3}\sqrt{\frac{a}{b}}$ there are two branches with different arrival times $t_1(v), t_2(v)$ both solving the same cubic. In this last case, only the smallest root contributes for $t_1\leq t \leq t_2$ and after $t_2$ both actually contribute.