Solving the following system of integrodifferential equations: speed up of the code

Consider a system of equations
$$
begin{cases}frac{partial f(E,t)}{partial t}-H[T_{gamma},f(E,t)]Efrac{partial f(E,t)}{partial E} -I[f(E,t),E,T_{gamma}(t)]=0, frac{partial T_{gamma}}{partial t} + frac{3H[T_{gamma},f(E,t)](4rho_{gamma}/3 + rho_{e}+ p_{e})+int limits_{0}^{infty} dE I[f(E,t),E,T_{gamma}(t)]}{partial rho_{gamma}(T_{gamma})/partial T_{gamma} + partial rho_{e}(T_{gamma})/partial T_{gamma}}=0end{cases}
$$
Here:
$$
I[f(E,t),E,T_{gamma}(t)] equivint limits_{Gamma(E)} dE’ G_{1}[E’,E,f(E,t),f(E’,t),T_{gamma}(t)]+ G_{2}(f(E,t),T_{gamma}(t),E)
$$
is some integral including $f(E)$, $f(E’)$ and $f(E)f(E’)$ terms, with $Gamma(E)$ being the domain of integration. $H,rho_{e},p_{e},rho_{gamma}$ are known functions; in particular,
$$
H = b sqrt{frac{pi^{2}}{30}g(T_{gamma}) + 6int limits_{0}^{infty}f(E)frac{E^{3}}{2pi^{2}}dE}, quad b = text{const},
$$
while $rho_{gamma},p_{e},rho_{e}$ do not include $f$.

Formally, $E$ varies in the limits $0<E<infty$, but below I restrict myself by the domain $0<Elesssim 30$.

The implementation of the solution of this system in Mathematica for a particular $I$, given below, follows an approach from this question:

(*Definitions*)
ME = 0.5;
etaBPlanck = 6.1*10^-10;
ng[T_] = 2*Zeta[3]/Pi^2 T^3;
mPl = 1.2*10^22;
hbar = 6.582119*10^-25;
sToGeVminusOne = 1/hbar;
sToMeVminusOne = sToGeVminusOne*10^-3;
rhoeTemp[T_] := 
 4/(2*Pi^2)NIntegrate[(Ee^2 Sqrt[Ee^2 - ME^2])/(
   Exp[Ee/T] + 1), {Ee, ME, Infinity}]
peTemp[T_] := 
 4/(6*Pi^2)NIntegrate[(Ee^2 - ME^2)^(3/2)/(Exp[Ee/T] + 1), {Ee, ME, Infinity}]
rhoe[T_] = 
 Quiet[10^Interpolation[
     Join[Table[{10^T, Log10[etaBPlanck*ng[10^T]*ME]}, {T, 
        Log10[10^-8], Log10[0.0049], 0.1}], 
      Table[{10^T, 
        Log10[etaBPlanck*ng[10^T]*ME + rhoeTemp[10^T]]}, {T, 
        Log10[0.005], Log10[19.99], 0.01}], 
      Table[{T, Log10[7/8*4*Pi^2/30*T^4]}, {T, 20, 300, 0.5}]], 
     InterpolationOrder -> 1][T]]
pe[T_] = Quiet[
  10^Interpolation[
     Join[Table[{10^T, -90}, {T, Log10[10^-8], Log10[0.0049], 0.1}], 
      Table[{10^T, Log10[10^-99 + peTemp[10^T]]}, {T, Log10[0.005], 
        Log10[19.99], 0.01}], 
      Table[{T, Log10[1/3 7/8*4*Pi^2/30*T^4]}, {T, 20, 300, 0.5}]], 
     InterpolationOrder -> 1][T]]
DrhoeDT[T_] = D[rhoe[T], T];
rhog[T_] = 2*Pi^2/30 T^4;
DrhogDT[T_] = D[rhog[T], T];


(*Integral I - continuous*)
dS11dE2[E1_, E2_] = -0.0034976546222801287` E1 E2^3*fn[E1] fn[E2];
S12[E1_, T_] = 0.02098592773368077` Exp[-E1/T] E1 T^4;
S21[E1_, T_] = -0.08394371093472308` E1 T^4 fn[E1];
dS221dE2[E1_, E2_, T_] = 
  0.03147889160052116`/(E1^2) Exp[-E1/
     T] T^2 (-E1^2 (E2^2 + 2 E2 T - 2 (-1 + Exp[E2/T]) T^2) - 
     2 E1 T (E2^2 + 2 E2 T - 2 (-1 + Exp[E2/T]) T^2) + 
     2 T^2 ((-1 + Exp[E2/T]) E2^2 - 2 (1 + Exp[E2/T]) E2 T + 
        4 (-1 + Exp[E2/T]) T^2)) fn[E2];
dS222dE2[E1_, E2_, T_] = 
  0.03147889160052116` /(E1^2) Exp[-E1/
     T] T^2 (2 (-1 + Exp[E1/T]) T^2 (E2^2 + 2 E2 T + 4 T^2) - 
     E1^2 (E2^2 + 2 E2 T - 2 (-1 + Exp[E1/T]) T^2) - 
     2 E1 T (E2^2 + 2 E2 T + 2 (1 + Exp[E1/T]) T^2)) fn[E2];
CollisionIntegralContinuous[E1_, T_] := 
 S12[E1, T] + S21[E1, T] + 
  NIntegrate[dS11dE2[E1, E2], {E2, 0, Infinity}] + 
  NIntegrate[dS221dE2[E1, E2, T], {E2, 0, E1}] + 
  NIntegrate[dS222dE2[E1, E2, T], {E2, E1, Infinity}]

(*Integral I - discretized*)
imax = 100;
DEvalue = 0.5;
Emin = 0.01;
EStep[k_, DE_] = Emin + DE*k;
fnThermalDiscrete[k_, Tn_] = Exp[-EStep[k, DEvalue]/Tn];
dS11dE2discrete[i_, j_] = 
  dS11dE2[E1, E2] /. {fn[E1] -> fn[i][t], 
     fn[E2] -> fn[j][t]} /. {E1 -> EStep[i, DEvalue], 
    E2 -> EStep[j, DEvalue]};
S12discrete[i_] = 
  S12[E1, T] /. {fn[E1] -> fn[i][t], fn[E2] -> fn[j][t]} /. {E1 -> 
     EStep[i, DEvalue], E2 -> EStep[j, DEvalue]};
S21discrete[i_] = 
  S21[E1, T] /. {fn[E1] -> fn[i][t], fn[E2] -> fn[j][t]} /. {E1 -> 
     EStep[i, DEvalue], E2 -> EStep[j, DEvalue]};
dS221dE2discrete[i_, j_] = 
  dS221dE2[E1, E2, T] /. {fn[E1] -> fn[i][t], 
     fn[E2] -> fn[j][t]} /. {E1 -> EStep[i, DEvalue], 
    E2 -> EStep[j, DEvalue]};
dS222dE2discrete[i_, j_] = 
  dS222dE2[E1, E2, T] /. {fn[E1] -> fn[i][t], 
     fn[E2] -> fn[j][t]} /. {E1 -> EStep[i, DEvalue], 
    E2 -> EStep[j, DEvalue]};
CollisionIntegrali[
  i_] := (S12discrete[i] + S21discrete[i] + 
    DEvalue*Sum[dS11dE2discrete[i, j], {j, 0, imax, 1}] + 
    DEvalue*Sum[dS221dE2discrete[i, j], {j, 0, i, 1}] + 
    If[i == imax, 0, 
     DEvalue*Sum[
       dS222dE2discrete[i, j], {j, i + 1, imax, 1}]]) /. {T -> Tg[t]}

(*A derivative dfn/dE - discretized*)
fnEnDerivative[i_] = 
  If[i != imax, (fn[i + 1][t] - fn[i][t])/DEvalue, -fn[imax][t]/
   DEvalue];
(*A function H*)
Hdynamical = 
  sToMeVminusOne/mPl Sqrt[
   8*Pi/3 (rhog[Tg[t]] + rhoe[Tg[t]] + 
      6*Sum[fn[i][t] 4*Pi (EStep[i, DEvalue])^3/(2*Pi)^3, {i, 0, imax,
          1}])];

(*System of equations and its solution*)
EquationsWithHubbleTable = 
  Join[{Tg'[t] + 
      1/(DrhoeDT[Tg[t]] + 
        DrhogDT[Tg[t]]) (4*Hdynamical*rhog[Tg[t]] + 
         3*Hdynamical*(rhoe[Tg[t]] + pe[Tg[t]]) + 
         DEvalue*Sum[
           CollisionIntegrali[i] 4*Pi EStep[i, DEvalue]^3/(
            2*Pi^3), {i, 0, imax, 1}]) == 0}, 
   Table[fn[i]'[t] - Hdynamical*EStep[i, DEvalue]*fnEnDerivative[i] - 
      CollisionIntegrali[i] == 0, {i, 0, imax, 1}]];
H[T_] = 1.66 Sqrt[10.75] T^2/mPl sToMeVminusOne;
t0 = 0.05;
Tt[t_] = Quiet[T /. Solve[H[T] == 1/(2 t), T][[2]]];
InitialConditionTable = 
  Join[{Tg[t0] == Tt[t0]}, 
   Table[fn[i][t0] == 1/(Exp[EStep[i, DEvalue]/(Tt[t0])] + 1), {i, 0, 
     imax, 1}]];
functionsTable = Join[{Tg}, Table[fn[i], {i, 0, imax, 1}]];
tmax = 4;
solWithDynamicalHubble = NDSolve[{EquationsWithHubbleTable, InitialConditionTable}, functionsTable, {t, t0, tmax}, Method -> {"EquationSimplification" -> "Solve"}][[1]];

My question is the following: is it possible to speed up the code above? In my machine, it requires 400 seconds or so to evaluate solWithDynamicalHubble for imax = 100 (apart from computing the table of equations), whereas my goal is to solve for imax = 20000.

P.S. My machine is 16 Gb RAM, i7-10875H. When running Mathematica, only ~15% of the CPU is used.