PT = {
{1/3, 1/3, 0, 1/3, 0, 0, 0, 0, 
            0}, 
{1/4, 1/4, 1/4, 0, 1/4, 0, 0, 0, 
            0}, 
{0, 1/3, 1/3, 0, 0, 1/3, 0, 0, 
            0}, 
{1/4, 0, 0, 1/4, 1/4, 0, 1/4, 0, 
            0}, 
{0, 1/5, 0, 1/5, 1/5, 1/5, 0, 1/5, 
            0}, 
{0, 0, 1/4, 0, 1/4, 1/4, 0, 0, 
            1/4}, 
{0, 0, 0, 1/3, 0, 0, 1/3, 1/3, 
            0}, 
{0, 0, 0, 0, 1/4, 0, 1/4, 1/4, 
            1/4}, 
{0, 0, 0, 0, 0, 1/3, 0, 1/3, 
            1/3}
};

P = Transpose[PT];

P // MatrixForm

x0 = {0, 0, 0, 0, 1, 0, 0, 0, 0};
x0 // MatrixForm

x1 = P . x0;
x1 // N // MatrixForm

x2 = P . x1;
x2 // N // MatrixForm

x3 = P . x2;
x3 // N // MatrixForm

Q = Array[q, 9];
Q // MatrixForm

sols = Solve[P . Q == Q]

Q = Q /. sols[[1]];
Q // MatrixForm

sumQ = Total[Q]

sols = Solve[sumQ == 1]

Q = Q /. sols[[1]];
Q // N // MatrixForm