Solution (source code)

= Solution

This is the <Blume–Capel model>. In the <mean-field approximation>, write $m=\langle\sigma_i\rangle$ and replace
$$
\sigma_i\sigma_j\simeq m\sigma_i+m\sigma_j-m^2.
$$
Because every site has <coordination number of a lattice>[coordination number] $q$, the resulting energy is
$$
E_{\rm MF}=\frac12NJqm^2+\sum_i\left[g\sigma_i^2-(Jqm+B)\sigma_i\right].
$$
The one-site <partition function> is therefore
$$
Z_1=\sum_{\sigma=-1}^{1}e^{-\beta[g\sigma^2-(Jqm+B)\sigma]}
=1+2e^{-\beta g}\cosh\!\left(\beta(Jqm+B)\right).
$$
With $\kappa=e^{-\beta g}$, the mean-field partition function is $Z_{\rm MF}=e^{-\beta NJqm^2/2}Z_1^N$. Taking $F=-T\log Z_{\rm MF}$ gives
$$
\boxed{\frac FN=\frac12Jqm^2-T\log\!\left[1+2\kappa\cosh\!\left(\beta(Jqm+B)\right)\right]}.
$$