How can I solve unintegrated Boltzmann equation fast and robustly?

Question

Consider the following integro-differential equations (the unintegrated Boltzmann equations for ultrarelativistic particles in the expanding Universe):

$$ \left(\frac{\partial}{\partial t} - H p \frac{\partial}{\partial p}\right) f(i,p,t) = C_i \left[f(i,P,t), p\right] \\ C_i = \int \mathrm{d}p_2 \mathrm{d}p_3 \, F(f(i,P,t),f(j,P,t)) M(p,p_2,p_3) \\ H(t) \sim \sqrt{\int \sum_i f(i,p,t)p^3 \, \mathrm{d}p} $$

These are the Boltzmann equations for some massless species i in an expanding Universe. There, f[i,p,t] is the distribution function of the species i; F[f[i,P,t],f[j,P,t]] is a polynomial in f[i], f[j] (with different momenta labels P), M is some polynomial in the momenta p,p2,p3, t is time.

How can I solve this system in Mathematica for a reasonable amount of time? Probably, it is some kind of discretization in p, but since I am a newbie, it would go extremely slow.

Explicit expressions for the system (using NIntegrate)

I have three species "ve", "vm", "vt".

Preliminary definitions for F (Fval[pt,pt1,pt,pt1,p1,p2,p3,p4,t]):

Do[Fval[pt, pt, pt, pt, p1_, p2_, p3_, p4_, 
   t_] = (1 - f[pt, p1, t])*(1 - f[pt, p2, t])*f[pt, p3, t]*
    f[pt, p4, t] - 
   f[pt, p1, t]*
    f[pt, p2, t]*(1 - f[pt, p3, t])*(1 - f[pt, p4, t]), {pt, {"ve", "vm", "vt"}}]
GF = 1.167*10^-5.;
Gval = 6.7*10^-39;
ffd[p_, T_] = 1/(Exp[p/T] + 1);
tT[T_] = 
 5.68*10^-25 (t /. 
    Refine[Solve[1/(2 t) == Sqrt[(8 Pi*Gval)/3 Pi^2/30 3*T^4], t], 
      T > 0][[1]])

The maximal momentum pmax to be considered, collision integrals Icoll[species1,species2,species3,species4,p,t] and the Hubble factor hubble[t]:

pmax = 1;
J1[p1_, p2_, p3_] = 
  16/15 p3^3 (10 (p1 + p2)^2 - 15 (p1 + p2)*p3 + 6*p3^2);
J2[p1_, p2_] = 16/15 p2^3*(10*p1^2 + 5*p1*p2 + p2^2);
J3[p1_, p2_, p3_] = 
  16/15 ((p1 + p2)^5 - 10 (p1 + p2)^2*p3^3 + 15*(p1 + p2)*p3^4 - 
     6*p3^5);
Table[With[{pt = pt, pt1 = pt1},
   Icoll[pt, pt1, pt, pt1, p1_, t_] := 
    GF^2/((2 Pi)^3 p1^2) (NIntegrate[
        Fval[pt, pt1, pt, pt1, p1, p2, p3, p1 + p2 - p3, t]*
         J1[p1, p2, p3], {p3, 0, p2}, {p2, 0, p1}] + 
       NIntegrate[
        Fval[pt, pt1, pt, pt1, p1, p2, p3, p1 + p2 - p3, t]*
         J2[p1, p2], {p3, p2, p1}, {p2, 0, p1}] + 
       NIntegrate[
        Fval[pt, pt1, pt, pt1, p1, p2, p3, p1 + p2 - p3, t]*
         J3[p1, p2, p3], {p3, p1, p1 + p2}, {p2, 0, p1}] + 
       NIntegrate[
        Fval[pt, pt1, pt, pt1, p1, p2, p3, p1 + p2 - p3, t]*
         J1[p1, p2, p3], {p2, p1, pmax},{p3, 0, p1}] + 
       NIntegrate[
        Fval[pt, pt1, pt, pt1, p1, p2, p3, p1 + p2 - p3, t]*
         J2[p2, p1], {p2, p1,pmax},{p3, p1, p2}] + 
       NIntegrate[
        Fval[pt, pt1, pt, pt1, p1, p2, p3, p1 + p2 - p3, t]*
         J3[p1, p2, p3], {p2, p1,pmax},{p3, p2, p1 + p2}])
   ], {pt1, {"ve", "vm", "vt"}}, {pt, {"ve", "vm", "vt"}}];
hubble[t_] := 
 Sqrt[(8*Pi*Gval)/3] Sqrt[
  NIntegrate[(f["ve", p, t] + f["vm", p, t] + f["vt", p, t]) p^3/(
    2 Pi^2), {p, 0, pmax}]]

Initial conditions and the final time for the system:

T0 = 3*10^-3;
t0 = tT[T0];
f0["vm", p_] = f0["vt", p_] = ffd[p, T0];
f0["ve", p_] = ffd[p, T0] + Exp[-(p - 0.05)^2/T0^2];
tmax = 10^3;