Partition Function for Rotational States

Do[
t = 
        10;    
z = (2  j + 1) 
          Exp[-j (j + 1)/t]; 
Print["j = ", 
        j, ", z= ", z // N], 
{j, 0, 10}]

Heat Capacity of Rotational States

Clear[t, j, z];   
jmax = 20;

Z = 
        Sum[(2  j + 1) Exp[-(j (j + 1)*b)], {j, 0, 
            jmax}];

energy = (-1/Z)*D[Z, b]  /.  b -> 1/t;

heatcap = D[energy, t];

plot1 = Plot[heatcap, {t, 0.01, 4}]