β=1/50. ;   (* try values from 1/1. (low temperature) to 1/30. (high temp) *)
ϵ = 1. ;       (* set rotation energy scale to 1 *)
ListPlot[
Table[
(2n+1) Exp[-n(n+1)β ϵ],
{ n,0,20 }
],
Filling -> Axis ,
PlotRange -> All
]

Clear[ β , ϵ ] ;   (* remove any existing values for these variables *)

nterms = 10 ;     (* number of terms retained in partition function sum *)

zz = Sum[ (2n+1)(E^(-n(n+1)β ϵ)) , {n,0,nterms} ] ;   (* symbolic partition function *)

(* note the substitution β -> 1/(kT), this makes avgE a function of T
  in preparation for differentating with respect to T. *)

avgE = - D[ Log[zz], β ] /. β -> 1/(k T) ;   (* average energy <E> *)

heatCapacity = D[ avgE, T ] ;  

 (* for plotting purposes, set k and ϵ to 1 in C(T) expression,
   this corresponds to plotting with a dimensionless temperature kT/ϵ = T *)
plotC = heatCapacity /. { k -> 1, ϵ -> 1 } ;

Plot[ 
plotC, 
{T,0,5}, 
PlotRange -> All,
Frame -> True,
FrameLabel -> { kT/ϵ, C/k}
]