Einstein Theory of Heat Capacity of Solid

cv[θ_,T_]:=3 R ((((θ/T)(((Exp[θ/(2 T)])/(Exp[θ/T]-1)))))^2)

R=8.314

Plot[cv[1000,T],{T,10,1000},PlotRange->{{0,1000},{0,25}}]

Manipulate[
Plot[cv[θ,T],{T,0,1000},AxesLabel->{"T(K)","Cv (J/mol)"},PlotLabel -> "Einstein Capacity",  PlotRange->{{0,1000},{0,25}}],
{θ,100,2500}]

For[i=0,i<4,i++,Print[i]]

fd[0]=0

For[i=1, i<41, i++, y=0.05 * i; fd[i]=N[3 (y^3)(∫_0^(1/y))(((x^4)(E^x))/((((E^x)-1))^2))\[DifferentialD]x];Print[y," ",1/y," ",fd[i]]]

For[i=1, i<41, i++, y=0.05 * i; fe[i]=((((1/y)(((Exp[1/(2 y)])/(Exp[1/y]-1)))))^2);Print[y," ",1/y," ",fe[i]]]

ListLinePlot[{Table[{0.05*i,3 R fd[i]}, {i,1,40,1}],Table[{0.05*i,3 R fe[i]}, {i,1,40,1}]},AxesLabel->{"T/θ","Cv (J/mol)"},PlotLegends->{"Debye","Einstein"}]

cvd[θ_,T_]:=3 R N[3 (((((T )/θ)))^3)(∫_0^(θ/T))(((x^4)(E^x))/((((E^x)-1))^2))\[DifferentialD]x]

cv[1000,100]

cvd[1000,100]

Plot[{cv[85,T],cvd[85,T]},{T,5,80},PlotPoints->15,MaxRecursion->0,PlotLegends->"Expressions",AxesLabel->{"T(K)","Cv (J/mol)"},PlotLabel -> "Einstein-Debye Capacity",  PlotRange->{{5,80},{0,25}}]

Eave[w_,T_]:=h w c * 100 /(Exp[h w c *100 /(k T)]-1)*na / 1000

h=6.62608*^-34

k=1.38065*^-23

c=2.99792458*^8

na = 6.022*^23

Manipulate[
Plot[Eave[w,T]/(na k T/1000),{T,0,2000},AxesLabel->{"T(K)","Eave / kT"},PlotLabel -> "Average Vibrational Energy",  PlotRange->{{0,2000},{0,1}}],
{w,10,3000}]