I’m trying to write the code of the first-order approximate solution by myself. I have a question: when using Shur decomposition, how to theoretically solve the covariance matrix of variables according to the optimal policy function matrix and state transition matrix.

You mean that you solved for the policy functions already and now want to compute theoretical moments? That involves solving the discrete Lyapunov equation. If you don’t care about numerical efficiency, there is an analytical formula involving Kronecker products. See Lyapunov equation - Wikipedia