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;