Excuse me each great god, and the following two matlab program, where we have a problem, running for a long time, not the results
Procedure 1
T=0:1:10 0;
Fplot (T, Ce)


The function C=Ce (T)

Syms T E
H=6.6261 * 10 ^ - 34; C=3 * 10 ^ 8; Me=9.10956 * 10 ^ - 31; Ne=8.5 * 10 ^ 28; KB=1.38065 * 10 ^ - 23.
A=c/800 f * 10 ^ 9;
W=2 * PI * f;

Ef=(h * c)/(8 * me * c ^ 2) * (3/PI) ^ (2/3) * ne ^ (2/3);
U=Ef * [1 - (PI) ^ 2/12 * (KB * T/Ef) ^ 2 + (PI) ^ 2/80 * (KB * T/Ef). ^ 4);

B=1/(KB * T);
Np=1/(exp (b. * (E - U)) + 1);
P=(8 * SQRT (2 * PI * me ^ (3/2) * SQRT (E))/(h ^ 3);
F1=np. * P. * E;
F2=np. * P;

C=ne * diff (Ep, T);
End

Program 2
T=0:1:10 0;
Fplot (T, F1)

The function F (T)=F1

H=6.6261 * 10 ^ - 34; C=3 * 10 ^ 8; Me=9.10956 * 10 ^ - 31; Ne=8.5 * 10 ^ 28; KB=1.38065 * 10 ^ - 23.
A=c/800 f * 10 ^ 9; E=1.6 * 10 ^ - 19; Z=1;

Syms T T
Ef=(h * c)/(8 * me * c ^ 2) * (3/PI) ^ (2/3) * ne ^ (2/3);
U=Ef * [1 - (PI) ^ 2/12 * (KB * T/Ef) ^ 2 + (PI) ^ 2/80 * (KB * T/Ef). ^ 4);
F=int (t ^ 0.5/(1 + exp (t - U/(KB * t))), t, 0, inf);
end