Much of the research which has been done on special kinds of matrix equations has almost the same structure. In this paper, an iterative approach is proposed which is more general and includes all such matrix equations. This approach is based on the global least squares (GL-LSQR) method. To do so, a linear matrix operator, its adjoint and a proper inner product are defined. Then we demonstrate how to employ the new approach for solving the general linear matrix equations system. When the matrix equations are consistent, a least-norm solution can be obtained. Moreover, the optimal approximate solution to a given group of matrices can be derived. Also, some theoretical properties and the error analysis of the new method are discussed. Finally, some numerical experiments are given to compare the new iterative method with some existent methods in terms of their numerical results and illustrate the efficiency of the new method.