  function x = Stationary(T)
% T is a circular stochastic matrix
% x is nx1 with nonnegative entries that sum to 1 and Tx = x.
[n,n] = size(T);
% Since T-I is column diagonally dominant, we will have P = I in PA = LU.
% So jsut do this...
[L,U] = lu(T - eye(n,n));
% U has this structure...,
%
%              x  x  0  0  0  x
%              0  x  x  0  0  x
%      U =     0  0  x  x  0  x
%              0  0  0  x  x  x
%              0  0  0  0  x  x
%              0  0  0  0  0  0

x = zeros(n,1);
x(n) = 1;
% Solve U(1:n-1:1:n-1)*x(1:n-1) + U(1:n-1,n)*1 = zeros(n-1,1) for x(1:n-1)... 
x(n-1) = -U(n-1,n)/U(n-1,n-1);
for k=n-2:-1:1
    x(k) = (-U(k,k+1)*x(k+1)-U(k,n))/U(k,k);
end
% Normalize...
x = x/sum(x);

% -2 if structure of U was not used
% -1 if structure of U was used but not in an efficient way
% -1 if x not normalized or is otherwise incorrect