Since $M^2$ is symmetric, it is diagonalizable with orthogonal eigenvectors. Also the eigenvalues are positive. Let $(v_i)$ be an orthonormal basis of eigenvectors for $M^2$ with eigenvalues $\lambda_i>0$. Now express $\mathbf 1$ in terms of the eigenvectors as $\sum b_iv_i$ and write $\mathbf x=\sum a_iv_i$. The equations reduce to $a_i(\lambda_i-K)=b_i$, where $K$ is the Rayleigh quotient $(\sum \lambda_i a_i^2)/(\sum a_i^2)$.

I would approach this by defining for each $t$, $\alpha_i(t)=b_i/(\lambda_i-t)$ and then $F(t)$ to be the Rayleigh quotient for the corresponding family of $(\alpha_i(t))$, that is $$ F(t)=\left(\sum \lambda_i \frac {b_i^2}{(\lambda_i-t)^2}\right)\Big / \left(\sum \frac{b_i^2}{(\lambda_i-t)^2}\right). $$ If for some $t$, $F(t)=t$, then you have a solution to your equations (namely $a_i=\alpha_i(t)$). Multiplying through by the denominators on both sides, you see that $F(t)$ is a continuous function of $t$, bounded between $\min\lambda_i$ and $\max\lambda_i$. Hence there is always a solution to the equation (by the intermediate value theorem).

But I don't think you should expect to be able to find an analytic solution, as this seems to be essentially the same as solving a typical degree $2n$ polynomial equation.

Since $M^2$ is symmetric, it is diagonalizable with orthogonal eigenvectors. Also the eigenvalues are positive. Let $(v_i)$ be an orthonormal basis of eigenvectors for $M^2$ with eigenvalues $\lambda_i\ge 0$. Now express $\mathbf 1$ in terms of the eigenvectors as $\sum b_iv_i$ and write $\mathbf x=\sum a_iv_i$. The equations reduce to $a_i(\lambda_i-K)=b_i$, where $K$ is the Rayleigh quotient $(\sum \lambda_i a_i^2)/(\sum a_i^2)$.

I would approach this by defining for each $t$, $\alpha_i(t)=b_i/(\lambda_i-t)$ and then $F(t)$ to be the Rayleigh quotient for the corresponding family of $(\alpha_i(t))$, that is $$ F(t)=\left(\sum \lambda_i \frac {b_i^2}{(\lambda_i-t)^2}\right)\Big / \left(\sum \frac{b_i^2}{(\lambda_i-t)^2}\right). $$ If for some $t$, $F(t)=t$, then you have a solution to your equations (namely $a_i=\alpha_i(t)$). Multiplying through by the denominators on both sides, you see that $F(t)$ is a continuous function of $t$, bounded between $\min\lambda_i$ and $\max\lambda_i$. Hence there is always a solution to the equation (by the intermediate value theorem).

But I don't think you should expect to be able to find an analytic solution, as this seems to be essentially the same as solving a typical degree $2n$ polynomial equation.