grid.m

function A = grid(n,m)
    if ~exist('m','var'), m=n; end
    A = kron(eye(n),path(m)) + kron(path(n),eye(m));
end

gridmodel.m

function J = gridmodel(n,m,edgeweight,nodeweight)
    J = grid(n,m) * edgeweight;
    J = J + eye(size(J,1)) * nodeweight;
    assert(all(eig(J) >= 0.))
end

gs_update.m

function G = gs_update(J)
    G = -tril(J) \ triu(J,1);
end

path.m

function A = path(n)
    A = diag(ones(1,n-1),1) + diag(ones(1,n-1),-1);
end

plot_spectral.m

function plot_spectral(A)
    [n,m] = size(A);
    L = diag(sum(A,2)) - A;
    [V,D] = eig(L);
    [y,i] = sort(diag(D));
    x = V(:,i(2));
    y = V(:,i(3));

    figure
    hold on
    for i=1:n
        for j=i+1:m
            if A(i,j), plot([x(i),x(j)],[y(i),y(j)],'b-'); end
        end
    end
    hold off
end