Abstract: Consider the numerical solution of coupled Sylvester matrix equations using Krylov subspace iterative methods are discussed in this paper. Such equations play a fundamental role in many systems and control applications. We present the preconditioned global full orthogonalization method and the preconditioned generalized minimal residual method for solving coupled Sylvester matrix equations. Some theoretical results are given. As numerical results show, it is essential to use preconditioning in association with Krylov subspace methods.