Discrepancy in Calculating Stationary Distribution of a Markov Chain

Question

I'm studying a Markov chain problem question 15 where a rat runs through a maze with the following transition matrix P:

p = np.array(
   [[0.  , 0.  , 1. , 0.  , 0.  , 0. ],
    [0.  , 0.  , 1. , 0.  , 0.  , 0. ],
    [0.25, 0.25, 0. , 0.25, 0.25, 0. ],
    [0.  , 0.  , 0.5, 0.  , 0.  , 0.5],
    [0.  , 0.  , 0.5, 0.  , 0.  , 0.5],
    [0.  , 0.  , 0. , 0.5 , 0.5 , 0. ]])

I'm trying to calculate the stationary distribution using the eigenvector method and the code below:

from scipy.linalg import eig 
transition_mat = p
S, U = eig(transition_mat.T)
stationary = np.array(U[:, np.where(np.abs(S - 1.) < 1e-8)[0][0]].flat)
stationary = stationary / np.sum(stationary)

The code gives me the stationary distribution as [1, 1, 4, 2, 2, 2]/12, which seems reasonable. However, when I run simulations using the code below, I get different results, such as [0.11111111, 0.11111111, 0.22222222, 0.22222222, 0.22222222, 0.11111111].:

import numpy as np
import random

def round():
    pos = random.randint(0, 5)
    for _ in range(1000):
        pos = np.random.choice(np.arange(0, 6), p=p[pos])
    return pos + 1

n = 800_000
res = []
for _ in range(n):
    res.append(round())

Additionally, when I use matrix exponentiation to compute the stationary distribution using the code below, I also get [0.11111111, 0.11111111, 0.22222222, 0.22222222, 0.22222222, 0.11111111]:

n = 1000000000
result = np.linalg.matrix_power(p, n)
z = np.dot(np.ones(3) / 3, result)

I would like to understand why these methods yield different results and whether the code I've provided is correctly calculating the stationary distribution. Any insights or suggestions on this issue would be greatly appreciated.