%inv. distr. of a stochastic matrix

%transition matrix
d=3;
P=[0.7 0.2 0.1;
   0.4 0.6 0;
   0   1   0];
%symbolic version
Ps=sym(P);

%powers
Q=eye(d); 
for (n=1:15)
    Q=Q*P;
    fprintf('P^%d=\n', n);
    for (i=1:d)
        for (j=1:d)
            fprintf('%6f ',Q(i,j));
        end
        fprintf('\n');
    end
    pause(0.5)
end

%left eigenvector for eigenvalue 1
  [v,e,vp]=eig(Ps');
 ind=(eval(diag(e))==1); 
 f1=find(ind==1); 
 f2=find(vp==f1); 
 v1=v(:,f2); 
 v1=v1/sum(v1);
 v2=eval(v1);

 fprintf('pi=[');
 for (i=1:d)
    fprintf('%s', rats(v2(i)));
 end 
 fprintf(']\n');
 

 
