Why is my theoretical answer of point of convergence and answer from simulations not the same?

Question

I am simulating a dynamic model which looks like the following:

$$ R(t+1) = AR(t)A' - \Gamma + I $$

The matrices $A, \Gamma, I$ are all 3x3 and known matrices. When I perform the simulation of this model, I notice that the values of $R$ converge. Therefore, I am trying to theoretically determine the equilibrium point. This means that I am trying to solve the following formula for $R^*$ (the equilibrium point of $R$):

$$ R^* = AR^*A' - \Gamma +I $$ For simplification $$ C := \Gamma +I, $$

Which gives $$ R^* = AR^*A' - C $$

I have tried to solve this formula by rewriting it in the form of a Sylvester equation which looks likes the following:

$$ ZX+XB = F $$

A closed-form solution for vec(X) in this equation is namely known. This would therefore allow me to determine the equilibrium point. I have rewritten the second equation (with the equilibrium point) in the following manner:

$$ B = A', X = R^*, Z = - A^{-1}, F = A^{-1}C $$

If you use these transformations in the Sylvester equation you achieve:

$$ - A^{-1}R^*+R^* A' = A^{-1}C $$ $$ A^{-1}R^* = R^*A' - A^{-1}C $$

If you then premultiply bij $A$ you get $$ R^* = AR^*A' -C $$

By being able to rewriting my equilibrium equation in a Sylvester equation, I was under the assumption that I could solve for $R^*$ by the solution to the Sylvester equation namely:

$$ (I\otimes Z+B'\otimes I)vec(X) = vec(F) $$ And then you can solve by

$$ vec(X) = (I\otimes Z+B'\otimes I)^{-1}vec(F), $$ by inserting the variables that I have defined earlier as B, X,Z and F. However, when I compare this theoretical result with my simulations result I get different answers for the convergence point. Therefore, I am wondering whether I made any mistake or miscalculation. Could you help me?

Edit: I have added the code to make the question more comprehensive. I have checked whether the answer that I get for R solves the equation ZX + XB = F and this is the case. So the mistake must be in the part where I transform the equilibrium formula into the Sylvester equation I think..

clc 
clear all
n=3;

A = [0 -0.4 0.4; -0.4 0 -0.4; 0.4 -0.4 0];
Gamma = [0 0 -1.6; 0 0 0; 0 0 -1.6];
I = eye(3);

R_etilde = zeros(n);

for t = 1:30
    R_etilde(:,:,t+1) = A*R_etilde(:,:,t)*A' - Gamma + I;
end

Z = -inv(A);
B = A'; 
F = inv(A)*(Gamma + I);
M = kron(I,Z)+kron(B',I);
vecF = F(:);    
vecR = inv(M)*vecF;
R = reshape(vecR, [3,3]);