Problem 6.21

Clear["Global`*"]

En[n_] := ϵ (1.03 n - 0.03 (n^2))

Start with just 2 terms to see what it looks like.

Z[β_] := Sum[Exp[-β  En[n]], {n, 1, 2}]

Z[β]

The average value of the energy is calculated symbolically by:

Ebar[β_] = ((-1)/(Z[β]))D[Z[β], β]

You can express that as a function of temperature T by using a replacement rule, replacing β by 1/(k T).

EbarT[T_] = Ebar[β] /. β -> 1/(k T)

Then the heat capacity is the derivative of the energy with respect to the temperature.

c[T_] = D[EbarT[T], T]

This is quite complicated, but can be simplified:

c[T_] = Simplify[D[EbarT[T], T]]

Calculate the heat capacity for a scaled temperature ‘t = k T / ϵ’, or T = ϵ t / k.
To get a dimensionless quantity, divide by ‘k’.

cscaled[t_] = c[ϵ t / k] / k

Plot[cscaled[t], {t, 0.1, 5}, PlotRange -> All]

Zideal[β_] := (1/(1 - Exp[-β ϵ]))

Eidealbar[β_] := ((-1)/(Zideal[β]))D[Zideal[β], β]

EidealbarT[T_] = Eidealbar[β] /. β -> (1/(k T))

cideal[T_] = Simplify[D[EidealbarT[T], T]]

cidealscaled[T_] = cideal[ϵ t / k] /k

Plot[{cscaled[t], cidealscaled[t]}, {t, 0.1, 5},PlotLegends -> {"nonideal", "ideal"},  PlotRange -> All, LabelStyle->Larger, AxesLabel->{"t", "c/k"}, ImageSize-> Scaled[0.7]]