NDSolve and NIntegrate with parameters

I have a problem with a numeric computation. I should solve a numeric ODE and the rhs is a numeric parametric function that I find with several NIntegrate. I’d like to fix the parameter (l, m is just a number) with a condition but when Mathematica tries to solve an integration gives me the error:
"the integrand U has evaluated to non-numerical values for all sampling points in the region with boundaries {0,500}".
The code is longer, I extrapolate the important things (I hope). Any help would be greatly appreciated!

m := 1000
theta := ArcSin[Sqrt[7*10^(-11)]]/2
Mp = 1.2209*10^22
L2[z_?NumericQ, x_?NumericQ] := Log[Abs[(((2 x)/z^2 + (2 Sqrt[m^2*z^2x^2])/z^2 + m^2) ((2 x)/z^2 - 2/z^2 Sqrt[m^2*z^2 + x^2] + m^2))/(((-2 x)/z^2 + (2 Sqrt[m^2*z^2 + x^2])/z^2 + m^2) ((-2 x)/z^2 - (2 Sqrt[m^2*z^2 + x^2])/z^2 + m^2))]]
L1[z_?NumericQ, x_?NumericQ] := Log[Abs[(-4 x/z^2 - m^2)/(4 x/z^2 - m^2)]]
Delta[z_?NumericQ] := (0.0071^2 z)/2
Sigma2[x_?NumericQ, z_?NumericQ, y_?NumericQ] := (l^4*x^3 (1 - y)^2)/((2  
[Pi])^3*4*z^5*(E^x + 1)*(((2*x*(1 - y))/z^2 - m^2)^2 + (l^4*m^4)/(32* 
[Pi])^2))
Sigmaint[x_?NumericQ, z_?NumericQ] := NIntegrate[Sigma2[x, z, y], {y, -1, 1}, 
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 1000},  MaxRecursion -> 20, AccuracyGoal -> 20]
Rate[z_?NumericQ] := NIntegrate[Sigmaint[x, z], {x, 0, 1000}, Method -> {"GlobalAdaptive"}, MaxRecursion -> 20, AccuracyGoal -> 20, Exclusions ->Automatic]
Gm[z_?NumericQ] := Rate[z]
Vtot[z_?NumericQ] := V[z]
G[z_] = 1/(1 + E^(z/z))
Gstar[z_] =If[z < 0.00001, 106.75, If[z < 10^(-4), 86.25, If[z < 10^(-3), 61.75, 10.75]]]
H[z_] = Sqrt[(4*Pi^3*Gstar[z])/(z^4*45*Mp^2)]
F[z_?NumericQ] := (Gm[z]*S2eff[z])/(4*H[z]*z)*G[z]
sol = ParametricNDSolve[{y'[z] == F[z], y[10^(-10)] == 0}, y, {z, 10^(-10), 10}, l]