In this paper, we develop a stationary iterative method to find the solution vector x for the invertible $n \times n$ linear system \[ (1.1)\qquad Ax = (I - B)x = f. \] (I and $I_k $ represent the appropriate $k \times k$ identity matrix.) We find x by replacing (1.1) with the equivalent system \[ (1.2)\qquad \tilde A\tilde x = (I - \tilde B)\tilde x = \tilde f,\quad \tilde x = x_0 \otimes x \] Solution vector x of (1.1) will be “easy to extract” from the solution of (1.2) since $\tilde x = x_0 \otimes x$ is always a decomposable tensor. For any quadratic polynomial where $\varphi (1) = 1$, we may construct tensor iteration matrix $\tilde B$ of (1.2) whose eigenvalues $ \pm \lambda $, are all determined by $\varphi $ according to the equation $( * )\lambda ^2 = \varphi (\mu )$ where $\mu $ runs over the eigenvalues of B in (1.1). With the ability to shape the spectrum of $\tilde B$ as per $( * )$, we develop an optimal stationary iterative algorithm to solve (1.2) in the special case when the spectrum of ...