clear all;

T=[0.4 0.3 0.2 0; 0.3 0.3 0.3 0.2; 0.2 0.2 0.3 0.5; 0.1 0.2 0.2 0.3];

x=[0; 1; 0; 0]; % state vector for day 0 (today)

X(:,1)=x;

for k=1:10
    X(:,k+1)=T*X(:,k);
end

X

[V, D]=eig(T);

Dlist=diag(D); % diag converts between a diagonal matrix and a column vector listing the diagonal elements

ind=find(abs(Dlist-1)<1e-7); % find the component index of the eigenvalue equal to 1 (within a certain tolerance)

Vscaled=V(:,ind)/sum(V(:,ind)) % scale the entries of the "ind"th column of V so that they become proportions/probabilities

