How to integrate a dependent function using symbolic integral methods (Shockley Queisser Limit)

I wish to take the integral of a function of form integral[((x^2)dx)/(e^x -1)]. The variable x here is term xg or xgc which is dependent on the variable Vg in the script below. The script I have included shows the working form of this code that I have ran correctly already using the trapz method (Intqs and Intqc), however in later calculations of this code the trapz method seems to not provide adequate/precise enough results. Ultimately I want to know what coding I need to input to integrate Qsfun and Qcfun listed below in the more symbolic form using integral(Qsfun, xg, inf) and integral(Qcfun, xgc, inf) with the limits from xg/xgc to inf respectively. I have tried it in several different ways and receive a multitude of error answers.

This is a portion of the Shockley Quiesser limit coding to determine the efficiency limit for single solar cell materials. My script is included here:

Vg = linspace(1.5e-1, 5, 3000);%Energy gap defined by us, allows optimization of where the cell has the highest efficiency.
for t = 1:length(Vg)

xg(t) = (q*Vg(t)) /(k*Ts); %Xg value operating on the energy gap value and dependent on the temp of the sun.
xgc(t) = (q*Vg(t)) /(k*Tc); %Xg value operating on the energy gap value and dependent on the temp of the cell.

x = linspace(xg(t),1000, 10000);%Size iteration of x based on values xg(t) to 1000, with 10000 interval points
for l = 1:length(x);

Qsfun(l) =(x(l)^2)/(exp(x(l))-1);
end
a = linspace(xgc(t),30*xgc(t), 30000);%Size iteration of a based on values xgc(t) to 30*xgc, with 30000 interval points
for y = 1:length(a);

Qcfun(y) =(a(y)^2)/(exp(a(y))-1);
end
Intqs(t) = trapz(x,Qsfun);
Intqc(t) = trapz(a,Qcfun);
%Iqs(t) = xg*trapz(x,Qsfun);
%Iqc(t) = xg*trapz(x,Qcfun);

%The number of quanta of frequency greater than vg incident per unit area time
Qs(t) = ((2*pi*((k*Ts)^3))/((h^3)*(c^2)))*Intqs(t);
Qc(t) = ((2*pi*((k*Tc)^3))/((h^3)*(c^2)))*Intqc(t); 

b = linspace(1e-10, 100, 3000);

for mm = 1:length(b)
Psfun(mm) = (b(mm)^3)/(exp(b(mm))-1); %Function of Plancks distribution integral
end

IntPs = trapz(b,Psfun);

u(t) = (xg(t)*Intqs(t))/IntPs;
%Plot of Vg vs u gives the dependence of the ultimate efficiency (u) upon the band gap(Vg) of the semiconductor.
%u(xg) = (xg*Intqs)/(IntPs); %Ultimate efficiency as a function of xg. Qs/Ps
%The Iqs function needs to be corrected by multiplying by xg according to eq. 2.8.