
P = [1/2 0 1/3; 1/4 3/4 1/3; 1/4 1/4 1/3];


pkg load symbolic   % skip this line on Matlab

syms lambda

B = P - lambda*eye(3)

det(B)
f = expand(det(B))

% by hand, note that 1 is an eigenvalue and divide
simplify(f / (lambda - 1))
% then find roots of the quadratic

% on computer, we can just
factor(f)

% same as:
eig(P)

% we determined the equilibrium prob vector
p = [2 6 3].'/11

% an experiment
x0 = rand(3, 1);  x0 = x0 / sum(x0);
x = x0;
for i=1:20
  x = P*x
end

% compare:
x
p
x - p   % getting close to zero
