# assign matrices
A <- matrix(c(.8, .4, .1, .5), 2, 2)
X <- matrix(c( 1, .5, .5,  2), 2)

# compute with above formulas
eg  <- eigen(A)
U   <- eg$vectors
U_t <- t(U)
las <- eg$values
T.  <- crossprod(U, X %*% U)
Z   <- T. / (1 - tcrossprod(las))
S   <- solve(U_t, t(solve(U_t, t(Z))))

# naive solution
nai <- X
for(i in 1:1000)
  nai <- crossprod(A, nai %*% A) + X

nai - S
#R           [,1]      [,2]
#R [1,] -3.55e-15 -4.44e-16
#R [2,] -8.88e-16  0.00e+00