代码：
h = 6.62607015*10^-34 ;(*普朗克常数*)
c = 299792458;(*光速*)
k = 1.380649*10^(-23);(*玻尔兹曼常数*)
e = 1.6*10^-19;
b[x_] := (1 - x) (11.4 - 19.4 x);
EG[x_] := x 0.64 e + (1 - x) 3.43 e - b[x] x (1 - x) e;
Subscript[EG, 1] = 3.43 e;
Subscript[EG, 2] = EG[0.5];
\[CapitalDelta]EG = (Subscript[EG, 1] - Subscript[EG, 2]);
(2 \[Pi] e^3)/(h^3 c^2)
NIntegrate[\[Epsilon]^2/(
E^((\[Epsilon] e)/(k 6000)) - 1), {\[Epsilon], Subscript[EG,
1], \[Infinity]}]
(2 \[Pi] e^3)/(h^3 c^2)
Integrate[\[Epsilon]^2/(
E^((\[Epsilon] e)/(k 6000)) - 1), {\[Epsilon], Subscript[EG,
1], \[Infinity]}]
另外这个积分里面怎么加一个待定的参数啊

这是因为普通的计算过程和数值的计算过程不同，普通的需要求出反导函数，所以简单函数可以用Intergrate,复杂的就不行。

应该算是bug，版本13.1返回的是未计算，并且这个积分的符号解其实可以算，只要把所有浮点数换成准确数：
test = Integrate[\[Epsilon]^2/(E^((\[Epsilon] e)/(k 6000)) - 1) //
Rationalize[#, 0] &, {\[Epsilon],
Subscript[EG, 1] // Rationalize[#, 0] &, \[Infinity]}]
N[test, 16]
(* 0.3336557962069662 *)