You have all the pieces... its just awkward to organize (below I use cells) and display.. Try something like
n = 100;
P = cell(1,n);
P{1} = [0.942 0 0 0 0.058 0 0 0; 0 0.982 0 0 0.018 0 0 0; 0 0 0.999 0 0 0.001 0 0; 0 0 0 0.993 0 0.007 0 0; 0 0 0 0 0.987 0 0.013 0; 0 0 0 0 0 0.999 0.001 0; 0 0 0 0 0 0 0.972 0.028; 0 0 0 0 0 0 0 1];
x_0 = [1000; 800; 300; 400; 700; 200; 300; 0];
% I expect x_0 to be a probability
x_0 = x_0/sum(x_0);
% Create P, P^2, P^3,...
for i = 2:n
P{i} = P{1}*P{i-1};
end
g = @(A)A*x_0;
evolvedX = cellfun(g,P,'un',0); % applies P^k to x_0 for k = 1,2,3,...
plot([evolvedX{:}]')
xlabel('Time-step')
ylabel('Probability')
legend('x_0(1)','x_0(2)','x_0(3)','x_0(4)','x_0(5)','x_0(6)','x_0(7)','x_0(8)')
