function [] = StationaryDist()
% Find the stationary distribution, two ways

p = [1/3, 1/3, 1/3 ; 7/10, 3/10, 0 ; 1, 0, 0];

v = null((transpose(p - eye(length(p))))); % the nullspace of p
dist1 = v.*(ones(length(p),1)*(1./sum(v))) % Changing normalization (null returns an orthonormal basis for p)


A = [-2/3, 1/3, 1 ; 7/10, -7/10, 1 ; 1, 0 ,1];

dist2 = [0,0,1]*inv(A)

end