function x = StationaryVec(A)
% A is an nxn matrix with positive entries that satisfy sum(A(:,k)) = 1
% for k=1:n.
% x is a vector with positive components that sum to one and satisfy Ax = x

[n,n] = size(A);
I = eye(n,n);
[L,U,P] = lu(A - I);
% A-I is singular so P(A-I) is singular so LU is singular so U is singular.
% This means thatU has a diagonal entry that is zero. It can be shown that
% U(n,n) is zero and all the other diagonal entries are nonzero.
% Find z(1:n-1) so U*[z;1] = 0. i.e., U(1:n-1,1:n-1)*z = -U(1:n-1,n)...
z = U(1:n-1,1:n-1)\(-U(1:n-1,n));
x = [z;1];
x = x/sum(x);