hi all! I've got something!
function u = proba(T)
k=1.380648813e-23;
u0=7.370012199e-19;
options=optimset('TolFun',1.0e-16,'TolX',1.0e-16);
u=fsolve(@(u)(2*u0^(3/2)/3-quad(@(e)(sqrt(e)./(exp((e-u)./(k.*T))+1)),0,1000000,1e-20)),u0,options);
plot(T,u);
hold on
end
only that...