clear all;

H=[.919 .175 0; .06 .825 .049; .021 0 .951];

x=[100000; 70000; 30000]; % state vector for year 0
X(:,1)=x;

for k=1:20
    %x=H*x;
    X(:,k+1)=H*X(:,k);
end

%format long g
%x
X(:,end-5:end)

[V, D]=eig(H);

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

Vscaled*sum(x) % multiply the proportions by the total number of households to get the number of households in each category