How to solve the following integrodifferential equation by generic approach?

Question

Consider the following equation: $$ \frac{\partial f}{\partial t} - p H(t)\frac{\partial f}{\partial p} = \mathcal{I}[p,t], \tag 1 $$ Here, f = f[p,t], with p being momentum and t time, while $H$ and $\mathcal{I}$ are

$$ H(t) = G \sqrt{\int p^{3}f dp}, $$ with $G$ being some constant, and $$ \mathcal{I}[p,t] = \sum_{i}\int \limits_{p_{2,\text{min}}^{i}}^{p_{2,\text{max}}^{i}}dp_{2}\int \limits_{p_{3,\text{min}}^{i}}^{p_{3,\text{max}}^{i}}dp_{3} \ \text{integrand}^{(i)} $$ where integrand depends on f with the various linear combination of momenta p, p2, p3.

In principle, the second summand in the LHS of this equation may be eliminated by the change of variable $p \to \epsilon = \frac{pa}{a_{0}T_{0}}$ ($a_{0} = 1, T_{0} = 0.015$ in units of GeV), where $a$ is defined as $$ H(t) \equiv \frac{1}{a}\frac{\partial a}{\partial t} $$ which would also allow the solution to stay within the same domain of the momentum coordinate $\epsilon$ at any moment in time: $$ \frac{\partial f(\epsilon,t)}{\partial t} = \mathcal{I}[\epsilon,t] $$ The system of equations represents a small piece of the equations from this question.

For the initial condition and boundaries, I take

t0 = 0.001*1/(6.58*10^-25);
pmin=0;
pmax=5;
fini[p_]=1/(Exp[p/0.015] + 1) + Exp[-30 (p - 0.1)^4];
fboundary[t_] = 0;

where fboundary is for p>pmax and p<0.

I wonder how it may be possible to solve this system of equations in Mathematica. The starting point for solving this system is probably the method of lines. However, there is the problem that the integral I is two-dimensional, and its discretization will produce hundreds of thousands of terms.

However, for sure, some more or less efficient way to solve the mentioned equation using discretization exists. In the literature, people use matlab, fortran, and python. In 2210.10307, for instance, they use matlab with ode15s solver, discretizing the collision integral by using Simpson's method (50-100 bins); the total time required to solve the differential equation is O(1 day), but this is for the coupled system of three equations of the type of $(1)$ with a much larger function $\mathcal{I}$, and also assuming a very good precision of the order of $10^{-4}$ in quantity $\int p^{3}f dp$ (I do not need it, maybe 10^-3 would be enough). 1512.02205 uses fortran, performing the momentum discretization. Neither of these codes is publicly available.

The integrands, the integration boundaries (6 in total), and the value of the constant G are provided below:

GF = 1.167*10^-5;
tabintegrands = {2/(15 Pi^3)
     GF^2 p1^(-2) p3^3 (10 (p1 + p2)^2 - 15 (p1 + p2) p3 + 
      6 p3^2) (-f[p1, t] f[p2, 
        t] (1 - f[p1 + p2 - p3, t]) (1 - f[p3, t]) + (1 - 
         f[p1, t]) (1 - f[p2, t]) f[p1 + p2 - p3, t] f[p3, t]), 
   2/(15 Pi^3)
     GF^2 p1^(-2) p2^3 (10 p1^2 + 5 p1 p2 + 
      p2^2) (-f[p1, t] f[p2, 
        t] (1 - f[p1 + p2 - p3, t]) (1 - f[p3, t]) + (1 - 
         f[p1, t]) (1 - f[p2, t]) f[p1 + p2 - p3, t] f[p3, t]), 
   2/(15 Pi^3)
     GF^2 p1^(-2) ((p1 + p2)^5 - 10 (p1 + p2)^2 p3^3 + 
      15 (p1 + p2) p3^4 - 
      6 p3^5) (-f[p1, t] f[p2, 
        t] (1 - f[p1 + p2 - p3, t]) (1 - f[p3, t]) + (1 - 
         f[p1, t]) (1 - f[p2, t]) f[p1 + p2 - p3, t] f[p3, t]), 
   2/(15 Pi^3)
     GF^2 p1^(-2) p3^3 (10 (p1 + p2)^2 - 15 (p1 + p2) p3 + 
      6 p3^2) (-f[p1, t] f[p2, 
        t] (1 - f[p1 + p2 - p3, t]) (1 - f[p3, t]) + (1 - 
         f[p1, t]) (1 - f[p2, t]) f[p1 + p2 - p3, t] f[p3, t]), 
   2/(15 Pi^3)
     GF^2 p1 (p1^2 + 5 p1 p2 + 
      10 p2^2) (-f[p1, t] f[p2, 
        t] (1 - f[p1 + p2 - p3, t]) (1 - f[p3, t]) + (1 - 
         f[p1, t]) (1 - f[p2, t]) f[p1 + p2 - p3, t] f[p3, t]), 
   2/(15 Pi^3)
     p1^(-2) ((p1 + p2)^5 - 10 (p1 + p2)^2 p3^3 + 15 (p1 + p2) p3^4 - 
      6 p3^5) (-f[p1, t] f[p2, 
        t] (1 - f[p1 + p2 - p3, t]) (1 - f[p3, t]) + (1 - 
         f[p1, t]) (1 - f[p2, t]) f[p1 + p2 - p3, t] f[p3, t])};
MapThread[(integrand[p1_, p2_, p3_, t_, #1] = #2) &, {{1, 2, 3, 4, 5, 
    6}, tabintegrands}];
pmax = 5;
tabboundaries[p1_, 
   p2_] = {{{0, p1}, {0, p2}}, {{0, p1}, {p2, p1}}, {{0, p1}, {p1, 
     p1 + p2}}, {{p1, pmax}, {0, p1}}, {{p1, pmax}, {p1, p2}}, {{p1, 
     pmax}, {p2, p1 + p2}}};
Gval = Sqrt[(3*6.7*10^-39)/(8*Pi)*(4*Pi)/(2*Pi)^3](**1/(6.58*10^-25)*);
Hintegrand[p_,t_]=Gval*Sqrt[p^3*f[p,t]];

The speed of evaluation of the integral I using fini is ~0.1 s (AdaptiveMonteCarlo method is used):

integral[p1_] := 
 Sum[NIntegrate[
   Evaluate[integrand[p1, p2, p3, t, i] /. f[P_, t] :> fini[P]], {p2, 
    tabboundaries[p1, p2][[i]][[1]][[1]], 
    tabboundaries[p1, p2][[i]][[1]][[2]]}, {p3, 
    tabboundaries[p1, p2][[i]][[2]][[1]], 
    tabboundaries[p1, p2][[i]][[2]][[2]]}, 
   Method -> "AdaptiveMonteCarlo"], {i, 1, 
   Length[tabboundaries[p1, p2]], 1}]
integral[RandomReal[{10^-5, pmax}]] // AbsoluteTiming
NIntegrate[Evaluate[Hintegrand[p,t]/.{f[p,t]->fini[p]}],{p,0,pmax}]