# How to solve a linear equation of qubit operators in qiskit?

I am trying to repeat the results from Peng and Kowalski (Eq 11) which solves the analytical form of an equation consisting of powers of the Hamiltonian (moments):

The $$K$$-th order PDS formalism is then associated with defining the introduced $$K$$ real parameters $$(a_1,\dots,a_K)$$ that minimize the value of $$R_K(s,a_1,\dots,a_K)$$. In this minimization process the necessary extreme conditions are given by the system of equations $$\frac{\partial R_K(s,a_1,\dots,a_K)}{\partial a_i} = 0, \hspace{4mm} (i=1,\dots,K),$$ which can be alternatively represented by the matrix system of equations for an auxiliary vector $$\textbf{X} = (X_1,\dots,X_K)^T$$ $$\textbf{MX}=-\textbf{Y}.$$ Here, the matrix elements of $$\textbf{M}$$ and vector $$\textbf{Y}$$ are defined as the expectation values of Hamiltonian powers (i.e. moments), $$M_{ij} = \langle \phi|H^{2K-i-j}|\phi\rangle,$$ $$Y_i = \langle \phi|H^{2K-i}|\phi\rangle$$ $$(i,j=1,\dots,K)$$.

Using the Jordan-Wigner transformation, I already get the Hamiltonian and moments in the qubit form but am having trouble finishing the solution. From above, $$X=(X_1,\dots, X_K)^T$$ ($$K=4$$ in my case) is the root waiting to be solved.