Not that it's directly relevant, but I have code for the generator matrix of a 1D Glauber-Ising model that could probably be reworked into 2D...

```
function y = glauber1d(symb,n,varargin);
% produces the generator matrix etc for a 1D Glauber-Ising model of n spins
% call as either glauber1d(1,n) for a less complete symbolic result, or
% glauber1d(0,5,[a,mu,H,kT,J]) for a more complete numerical result--i.e.,
% symb is a flag indicating whether or not to use symbolic calculations
% (this requires the symbolic toolbox in order to work)
% a (Glauber's alpha) is the spin flip rate, depends on the coupling
% between the GI system and the bath;
% mu is the magnetic moment associated with the spins;
% H is the magnetic field strength;
% kT is (well, you know);
% J is the exchange energy
if symb % SYMBOLICS
syms a b g real;
else % NUMERICS
args = varargin{1};
a = args(1);
mu = args(2);
H = args(3);
kT = args(4);
J = args(5);
b = tanh(mu*H/kT); % Glauber's beta (NOT 1/kT)
g = tanh(2*J/kT); % Glauber's gamma
end
% produce an array with rows equal to spin configurations
temp = dec2bin(0:((2^n)-1),n);
for j = 1:2^n
for k = 1:n
s(j,k) = 2*str2num(temp(j,k))-1;
end
end
% obtain spin flip rates
for j = 1:2^n
for k = 1:n
km = mod(k-2,n)+1;
kp = mod(k,n)+1;
temp = (g/2) * (b - s(j,k)) * (s(j,km) + s(j,kp));
w(j,k) = (a/2) * (1 - b*s(j,k) + temp);
end
end
% generator matrix
if symb
Q = sym(zeros(2^n));
else
Q = zeros(2^n);
end
for j1 = 1:2^n
for j2 = 1:2^n
if sum(abs( s(j1,:) - s(j2,:) )) == 2 % single spin flip
% now find out which spin gets flipped
k0 = find( s(j1,:) - s(j2,:) );
Q(j1,j2) = w(j1,k0);
end
end
end
if symb
Q = simplify( Q - diag(sum(Q,2)) );
else
Q = Q - diag(sum(Q,2));
end
% invariant distribution p (if you want it)
if 2^n - 1 - rank(Q)
'error'
y = 0;
return;
else
p0 = null(Q')';
end
if symb, simplify(p0); end
sp0 = sum(p0);
if symb, simplify(sp0); end
p = p0 / sp0; % invariant distribution
y = Q;
```