In this paper, a new iterative solution method is proposed for solving multiple linear systems. The proposed method is based on the global least squares (GL-LSQR) method. A linear operator is defined to connect all the linear systems together. To approximate all numerical solutions of the multiple linear systems simultaneously, the GL-LSQR method is applied for a new operator and the approximate solutions are obtained recursively.
The presented method is compared with the well-known LSQR method. Finally, numerical experiments on test matrices
are presented to show the efficiency of the new method.