(*this Mathematica notebook evaluates several observables of the 1D Ising model and plots them 
in the temperature/field plane, derivatives (of the free energy) are evaluated using the Mathematica functions 
D[f[x],x] for df/dx and Derivative[0,2][f[x,y]] for d^2f/dy^2 *)

(*1D Ising model with J=1,hence Jtilde=1/T=b *)

(*the two eigenvalues of the transfer matrix as a function of inverse temperature b and magnetic field*)
lambda1[b_,h_]:=Exp[b] Cosh[b h] + Exp[b] Sqrt[Sinh[b h]^2 + Exp[-4 b] ]
lambda2[b_,h_]:=Exp[b] Cosh[b h] - Exp[b] Sqrt[Sinh[b h]^2 + Exp[-4 b] ]

(*(ln Z)/N is equal the logarithm of the maximum eigenvalue, defined here explicitly*)
loglambda1[b_,h_]:=Log[lambda1[b,h]]

(*magnetisation evaluated by taking the derivative of lnZ wrt to h*)
m[b_,h_]=1/b D[Log[lambda1[b,h]],h]

(*magnetisation plotted in the T,h plane*)

Plot3D[m[1/T,h],{T,0,5},{h,-2,2}]

(* nearest-neighbour correlations in t,h plane*)
corr[b_,h_]:=Derivative[1,0][loglambda1][b,h]

Plot3D[corr[1/T,h],{T,0,5},{h,-2,2}]

(* zero-field nearest-neighbour correlations against temperature *)
Plot[corr[1/T,0],{T,0,5}]

(* zero-field nearest-neighbour correlations against inverse temperature *)
Plot[corr[b,0],{b,-.0001,5}]

(* - beta f = Log[lambda1]. plot f free energy in t,h plane*)

Plot3D[-T Log[lambda1[1/T,h]],{T,0,5},{h,-2,2}]

(* s=ln Z +beta e. entropy in t,h plane*)
(*the entropy per spin is ln2 at high temperatures, but reaches zero at low temperatures*)

Plot3D[Log[lambda1[1/T,h]]+1/T energy[1/T,h],{T,0,5},{h,-2,2}]

(*the magnetic susceptibility dm/dh is second derivative of lnZ wrt h*)
susceptibility[b_,h_]:= b Derivative[0,2][loglambda1][b,h]

(*note the divergence of the magnetic susceptibility at h=0 as the temperature is lowered*)
Plot3D[susceptibility[1/T,h] ,{T,0,5},{h,-2,2}]

(*the magnetic susceptibility at zero field*)
Plot[susceptibility[1/T,0] ,{T,0.,5}]

(*the specific heat is second derivative of lnZ wrt beta*)
specificheat[b_,h_]:=b^2 Derivative[2,0][loglambda1][b,h]

Plot3D[specificheat[1/T,h] ,{T,0,5},{h,-2,2}]

(*susceptibility at zero field*)

Plot[specificheat[1/T,0] ,{T,0,5}]

(*the difference between the two eigenvalues vanishes with decreasing temperature*)

Plot[{lambda1[1/T,0],lambda2[1/T,0]},{T,0,5}]