Z = Exp[-β E1] + Exp[-β E2]

Ebar = (E1 Exp[-β E1] + E2 Exp[-β E2])/Z

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

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