function pi = func_stationaryLaw(Q)
% Computes the stationary law of the transition matrix Q.
% Input
%   Q  k*k  transition matrix
% Output
%   pi 1*k  stationary law

    k = size(Q,1);
    
    QA = [(Q - eye(k)) ones(k,1)];
    sol = zeros(k+1,1);
    sol(k+1) = 1;
    pi = ((QA')\sol)';
