Heat capaity of a non-interacting J=1 system in a magnetic field.

Δ  =  g(Subscript[μ, B]) 
        H and introducing the reduced temperature
t = (((Subscript[k, B]) T)/Δ) = (1/β)Δ,  
    and the reduced heat capacity
\(c = \(C\_V\/\(\(k\_B\) N\)\ \ then\ \ c = \(c[
            t] = \(1\/t\^2\) \(e\^\(1/t\) + e\^\(\(-1\)/t\) + 4\)\/\((1 + e\^\(1/t\) + e\^\(\(-1\)/t\))\)\^2\)\)\)

c[
          t_] = (1/(t^2)) (((E^(1/t)) + (E^(-1/t)) + 4)/(((1 + (E^(1/t)) + (E^(-1/t))))^2));

g1 = 
      Plot[c[t], {t, 0, 7}, 
        AxesLabel -> {"t", "Heat capcity (c)"}, 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 14}, PlotStyle -> {Thickness[0.012]}]

clow[t_] = (1/(t^2)) (E^(-1/t));  
    chigh[t_] = (2/(3  (t^2)));
g2 = 
      Plot[clow[t], {t, 0, 0.6}, PlotStyle -> {Thickness[0.008]}, 
        DisplayFunction -> Identity]; 
    g3 = Plot[chigh[t], {t, 1.2, 7}, PlotStyle -> {Thickness[0.008]}, 
        DisplayFunction -> Identity];
Show[g1, g2, g3, DisplayFunction -> $DisplayFunction]

c[t]->(1/(t^2))(e^(-1/t)) in the limit of t->0, which expression vanishes exponentially. 
The high-temperature result is c[t] ->(2/(3  (t^2))).
The position of the maximum is:

FindRoot[c'[t] == 0, {t, 0.5}]

(Subscript[t, 0]) = t  /. %; c[(Subscript[t, 0])]

\(This\ behaviour\ of\ the\ heat\ capacity\ is\ called\ a\ Schottky\ anomaly . \ \ The\ number\ of\ levels\ involved\ in\ the\ anomaly\ is\ determined\ by\ \ the\ integral\ of\ C/
        T\ from\ zero\ to\ \(infinity : 
The\ reduced\ entropy\ \ is\)\)
s  = s[
          t] = (S/((Subscript[k, B]) N)) = (∫_0^t) c[
                      u]/u du  ->  Log[
                  2  J + 1] for t -> ∞
\(and\ may\ be\ calculated\ directly\ , \ \(\(but\)\(\ \)\(Mathematica\)\(\ \)\(runs\)\(\ \)\(into\)\(\ \)\(troubles\)\(\ \)\(if\)\(\ \)\(you\)\(\ \)\(try\)\(\ \)\(to\)\(\ \)\(do\)\(\ \)\(that . \ 
Here\)\(\ \)\(it\)\(\ \)\(is\)\(\ \)\(determined\)\(\ \)\(from\)\(\ \)\(the\)\(\ \)\(partition\)\(\ \)\(function\)\(\ \)\(as\)\(\ \)\(follows\)\(\ \)\)\)

z[t_] = 1 + (E^(1/t)) + (E^(-1/t));  
    s[t_] = Log[z[t]] + ((t z'[t])/z)[t];
Plot[s[t], {t, 0, 4}, AxesLabel -> {"t", "Entropy (s)"}, 
      TextStyle -> {FontFamily -> "Times", 
          FontSize -> 14}, PlotStyle -> {Thickness[0.012]}]

Limit[s[t], t -> ∞]

% // N

The result is Log[2J+1], where 2J +1 (here 3) is the number of levels.

Clear[c, s, t, z]