In this paper, a modified conjugate gradient iterative algorithm is proposed to solve the coupled Sylvester matrix equations \(\sum \nolimits _{j=1}^{q}A_{ij}X_{j}B_{ij}, (i\in \overline{1,p})\) . A novel theoretical proof is provided to confirm that the iterative sequence generated by the proposed algorithm converges to the solution of the considered matrix equations. Furthermore, it is proved that the proposed algorithm can generate some orthogonal matrix groups, allowing the solution of the matrix equations to be found within finite iteration steps. Additionally, the unique least Frobenius norm solution of the matrix equations is derived by choosing some special initial values. Finally, some numerical simulations are taken to substantiate the convergence of the proposed algorithm.
