(*this Mathematica notebook evaluates the free energy function of the Weiss ferromagnet (mean field ferromagnet) with J=1 *)

(* Z = exp{-N beta f(beta,h)} = \int dm exp{-N beta f(beta,h,m)}}    *)
(* we consider beta f(beta,h,m), denoted by betaf*)

betaf[beta_,h_,m_]:=beta m^2/2 - Log[2 Cosh[beta m + beta h]]

(* first we plot betaf as a function of m in the high temperature phase, 
beta=0.9<beta_c=1 *)

Plot[betaf[.9,0,m],{m,-1,1}]

(*at the critical temperature beta=beta_c=1*)
Plot[betaf[1,0,m],{m,-1,1}]

(*in the low temperature phase, beta=1.1>beta_c=1*)
Plot[betaf[1.1,0,m],{m,-1,1}]

(*explore the free energy function yourself by tuning beta and h.
Initially beta is set to the critical value beta=1 and h=0*)
(*for instance, what happens if you turn on a magnetic field, either below T_c or above T_c ? *)

Animate[Plot[betaf[beta,h,m],{m,-1,1}],{{beta,1},0.1,2},{{h,0},-.5,.5
},AnimationRunning->False]