function Aeq = is_stat(m,n)
% is the distribution P stationary?
mn = m^n;

Aeq = zeros(0,mn);
rind = 0;
for i=1:n-1
  j = i+1;
  % enforce constraint: marg(X_i) = marg(X_j)
  for y = 1:m
    rind = rind + 1;
    %Aeq(rind,1) = 0;
    for xi = 1:m^n
      x = i2c(xi,m,n);
      if x(i) ~= x(j)
        if (x(i)==y)
          Aeq(rind,xi) = 1;
        elseif (x(j)==y)
          Aeq(rind,xi) = -1;
        end
      end
    end
  end
end

return


function c=i2c(i,m,n)
dims = repmat(m,[1 n]);
c = ind2coord(i,dims);
return