mpl = 1/Sqrt[6.70837*10^-39];
gsT = 106.75;
Sup[\[CapitalLambda]I_?NumericQ, \[CapitalGamma]I_?NumericQ] :=
 Block[{a, \[Rho]r, Trad, tf, t},
  tf = 10/\[CapitalGamma]I;
  {a, \[Rho]r} =
   NDSolveValue[{a'[t] ==
      a[t]*Sqrt[(8 \[Pi])/(3 mpl^2) (\[Rho]r[
            t] + \[CapitalLambda]I^4/
             a[t]^3 Exp[-\[CapitalGamma]I t])], \[Rho]r'[t] +
       4*Sqrt[(8 \[Pi])/(3 mpl^2) (\[Rho]r[
             t] + \[CapitalLambda]I^4/
              a[t]^3 Exp[-\[CapitalGamma]I t])] \[Rho]r[
         t] == \[CapitalGamma]I \[CapitalLambda]I^4/
        a[t]^3 Exp[-\[CapitalGamma]I t],
     a[0] == 1, \[Rho]r[0] == 0}, {a, \[Rho]r}, {t, 0, tf},
    MaxStepFraction -> 10^-5, MaxSteps -> 10^6];
  Trad[t_] := (30/(\[Pi]^2 gsT) \[Rho]r[t])^(1/4);
  Trad[tf]*a[tf]
  ]