Quit[];

Para calcular funciones de partición de la red cuadrada. [Beale,  Exact Distribution of Energies in the Two-Dimensional Ising Model, PRL 76 78]

n = 8;  (*dimensión de la red. Para una red de 20x20 puede tardar unos 10 segundos, dependiendo de la computadora. Si van a trabajar con redes muy grandes, convendría que guardaran los resultados en un archivo para no tener que calcularlos cada vez*)

m=n; (*m = n para la red cuadrada. También pueden calcularse las funciones de partición de redes rectangulares*)

Clear[a,b,c,s,c2,s2,x,z,coef,z1,z2,z3,z4];
acc=Floor[N[n m Log[2]/Log[10]]( 3/2)];
If[Mod[n,2]==0,z1=2^(m n/2-1) Product[c2[k],{k,1,n-1,2}];
z2=2^(m n/2-1) Product[s2[k],{k,1,n-1,2}];
z3=2^(m n/2-1) c[0] c[n] Product[c2[k],{k,2,n-2,2}];
z4=2^(m n/2-1) s[0] s[n] Product[s2[k],{k,2,n-2,2}],z1=2^(m n/2-1) c[n] Product[c2[k],{k,1,n-1,2}];
z2=2^(m n/2-1) s[n] Product[s2[k],{k,1,n-1,2}];
z3=2^(m n/2-1) c[0] Product[c2[k],{k,2,n-1,2}];
z4=2^(m n/2-1) s[0] Product[s2[k],{k,2,n-1,2}]];
z1+z2+z3+z4;
c[0]=((1-x)^m+(x(1+x))^m)/2^(m/2);
s[0]=((1-x)^m-(x(1+x))^m)/2^(m/2);
c[n]=((1+x)^m+(x(1-x))^m)/2^(m/2);
s[n]=((1+x)^m-(x(1-x))^m)/2^(m/2);
b=2x(1-x^2);
a[k_]=(1+x^2)^2-b Cos[Pi k/n];
c2[k_]=2^(1-2m)(Sum[m!/j!/(m-j)!(a[k]^2-b^2)^(j/2) a[k]^(m-j),{j,0,m,2}]+b^m);
s2[k_]=2^(1-2m)(Sum[m!/j!/(m-j)!(a[k]^2-b^2)^(j/2) a[k]^(m-j),{j,0,m,2}]-b^m);
z=Rationalize[Chop[Expand[N[z1+z2+z3+z4,acc]]]];
coef=Table[0,{m n+1}];
Do[coef[[j+1]]=Coefficient[z,x,2j],{j,0,m n}];

(*el resultado se guarda en la lista coef*)

coef

Total[coef]
(2^(n^2))

La función de partición, la energía y el calor específico:

Z[x_]:= Sum[coef[[i+1]] (x^((2(n^2)-4i))), {i, 0, Length[coef]-1}];
U[x_]:= (1/((n^2)Z[x])) Sum[(2(n^2)-4i) coef[[i+1]] (x^((2(n^2)-4i))), {i, 0, Length[coef]-1}];
Cv[x_]:=((Log[x])^2) (  (1/((n^2)Z[x])) Sum[(((2(n^2)-4i))^2) coef[[i+1]]  (x^((2(n^2)-4i))), {i, 0, Length[coef]-1}]-(n^2) ((U[x])^2));

ListPlot[Table[{T, Cv[(E^(1/T))]}, {T, 0.1, 5., 0.05}]]

ListPlot[Table[{T, U[(E^(1/T))]}, {T, 0.1, 5., 0.05}]]