MathOverflow will be down for maintenance for approximately 3 hours, starting Monday evening (06/24/2013) at approximately 9:00 PM Eastern time (UTC-4).

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;