# define transition matrix
P = t(matrix(c(0.8, 0.2, 0.0, 0.4, 0.4, 0.2, 0.2, 0.6, 0.2), nrow = 3, ncol = 3))

# eigenvalues and eigenvectors
eig = eigen(P)

# get eigenvalues
v = eig$values

# diagonal matrix with eigenvalues
diagonal = diag(v)

# get eigenvectors
phi = eig$vectors

# load MASS package for matrix inversion
library(MASS)
phiINV = ginv(phi)
# perform matrix multiplication
phihoch176 = phi %*% diagonal^176 %*% phiINV

# multiply with initial distribution alpha = c(1,0,0)
c(1,0,0) %*% phihoch176