function pv=pvlu(P)
% Calcola il vettore invariante di probabilita' della matrice stocastica P
% mediante fattorizzazione LU di I-P
n=size(P,1);
A=eye(n)-P;
% Eliminazione gaussiana
for k=1:n-1
% Aggiorna la matrice L
    A(k+1:n,k)=A(k+1:n,k)/A(k,k);
% Aggiorna la matrice U
    for i=k+1:n
        A(i,k+1:n)=A(i,k+1:n)-A(i,k)*A(k,k+1:n);
    end
end
% Back substitution
pv=ones(n,1);
for j=n-1:-1:1
    s=0;
    for i=j+1:n
        s=s+pv(i)*A(i,j);
    end
    pv(j)=-s;
end
pv=pv/sum(pv);

    