#10 might be caused due to the distribution function f taking on negative values. One possible way to overcome this is to evolve log(f), which has no bounds instead of f, which is positive.
The equation for F(x, v, t) == log(f) has the same structure as f:
dF/dt = \partial F/\partial t + \partial F / \partial * dx/dt + \partial F / \partial v * dv/dt
where the characteristic dv/dt now depends on f == exp(F)