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

Question

Consider a system of equations $$ \begin{cases}\frac{\partial f(E,t)}{\partial t}-H[T_{\gamma},f(E,t)]E\frac{\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)](4\rho_{\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}}=0\end{cases} $$ Here: $$ I[f(E,t),E,T_{\gamma}(t)] \equiv\int \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}) + 6\int \limits_{0}^{\infty}f(E)\frac{E^{3}}{2\pi^{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<E\lesssim 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.0034976546222801287E1 E2^3*fn[E1] fn[E2]; S12[E1_, T_] = 0.02098592773368077 Exp[-E1/T] E1 T^4;
S21[E1_, T_] = -0.08394371093472308E1 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 + DEk;
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] + 
    DEvalueSum[dS11dE2discrete[i, j], {j, 0, imax, 1}] + 
    DEvalueSum[dS221dE2discrete[i, j], {j, 0, i, 1}] + 
    If[i == imax, 0, 
     DEvalueSum[
       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[
   8Pi/3 (rhog[Tg[t]] + rhoe[Tg[t]] + 
      6Sum[fn[i][t] 4Pi (EStep[i, DEvalue])^3/(2Pi)^3, {i, 0, imax,
          1}])];
(System of equations and its solution)
EquationsWithHubbleTable = 
  Join[{Tg'[t] + 
      1/(DrhoeDT[Tg[t]] + 
        DrhogDT[Tg[t]]) (4Hdynamicalrhog[Tg[t]] + 
         3Hdynamical(rhoe[Tg[t]] + pe[Tg[t]]) + 
         DEvalueSum[
           CollisionIntegrali[i] 4Pi EStep[i, DEvalue]^3/(
            2Pi^3), {i, 0, imax, 1}]) == 0}, 
   Table[fn[i]'[t] - HdynamicalEStep[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.