syms p T nH alpha2 density_part_integral density_photons_T density_photons_total

p=(1-nH/density_photons_total)*(density_photons_T/density_photons_total)

h=6.626*10^(-34);
kb=1.3806*10^(-23);
c=3*10^8;
V=10^(-6);

alpha2=8*pi*V*(kb*T/(h*c))^3;
fdensityT=x^2/(exp(x)-1);
fnumT=matlabFunction(fdensityT);
density_part_integral=integral(fnumT,0,13.6);
density_photons_T=alpha*density_part_integral;
density_photons_total=alpha*2*zeta(3);

solve(p==(1-nH/density_photons_total)*(density_photons_T/density_photons_total),T);