function integrand, x
   f = (x * x * x * x * exp(x))/((exp(x)-1.)*(exp(x)-1.))
   return, f
end

k = 1.381E-23
debyeT = 309.
T = 350.
xT = 0.
yS = 0.
while(T le 1083.) do begin
   uLimit = debyeT / T
   integral = qsimp('integrand', 0., uLimit, /DOUBLE,EPS=0.001)
   S = (9. * k * T * T * T / (debyeT * debyeT * debyeT)) * integral
   xT = [xT, T]
   yS = [yS, S]
   T = T + 100.
Endwhile

plot,xT,yS

end