In this research, we extend three attractive iterative methods—conjugate gradient, conjugate residual, and minimal residual—to solve large sparse symmetric multilinear system \(\mathcal {A}\textbf{x}^{m-1}=\textbf{b}\) . We prove that the developed iterative methods converge under some appropriate conditions. As an application, we applied the proposed methods for solving the Klein–Gordon equation with Dirichlet boundary condition. Also, comparing these iterative methods to some new preconditioned splitting methods shows that, applying new methods for solving symmetric tensor equation \(\mathcal {A}\textbf{x}^{m-1}=\textbf{b}\) , in which the coefficient tensor is an \(\mathcal {M}\) -tensor, are more efficient. Numerical results demonstrate that our methods are feasible and effective for solving this type of tensor equations. Finally, some concluding remarks are given.
