Explicit solution to a Rayleigh quotient equation For 5 months! I have been struggling to solve the following equations analytically without numeric method (ie, Newton method):

Main equation:
$$
 \biggl(M^2-\cfrac{\mathbf{x^{\text{T}}}M^2\mathbf{x}}{\mathbf{x^{\text{T}}}\mathbf{x}}E\biggr)\mathbf{x}=\mathbf{1}
$$
Constraint equations:
$$
\begin{cases}
 \mathbf{x^{\text{T}}1}=0 \\
\\
\mathbf{x^{\text{T}}x}=u 
\end{cases} $$
where $\{M,E\}\in\mathbf{R}^{n \times n}$ and $\{\mathbf{1},\mathbf{x}\}\in\mathbf{R}^n$ are defined, then $M$ is an arbitrary symmetric matrix, $E$
  is an identical matrix,
  $\mathbf{1}$ is all one vector, $\mathbf{x}$ is a
  variable vector and $u\in\mathbf{R}$ is a scalar.
  Furthermore, as a knowledge, the below equation form is called Rayleigh
  quotient $R(M^2,\mathbf{x})$:
$$R(M^2,\mathbf{x}):=\cfrac{\mathbf{x^{\text{T}}}M^2\mathbf{x}}{\mathbf{x^{\text{T}}}\mathbf{x}}$$

Now, we attempt to estimate the $\mathbf{x}$. Does the analytic solution or method exist? My ability is shortage but, I guess that this problem has a beautiful solution. Also, main equation is a simultaneous cubic equation. Theoretically, this is solvable. Just, this is my theme question.
 A: You may already know this, but your Main Equation is equivalent to the system
$$ \begin{cases} (M^2-\lambda E)\mathbf{x} = \mathbf{1} & (*) \\ \mathbf{x}^T (M^2 - \lambda E) \mathbf{x} = 0 & (**) \end{cases}. $$
The first constraint, $\mathbf{x}^T \mathbf{1} = 0$, follows from the equations, while the second one, $\mathbf{x}^T \mathbf{x} = u$. in general does not. An obvious remark then is that your equation may not have any solutions (unless the matrix $M$ is somehow special), because the number of equations exceeds the number of unknowns.
Here's a way to simplify the above system, under the assumption that $\lambda$ is not an eigenvalue of $M^2$. Let $W(\lambda) = (M^2-\lambda E)^{-1}$, which by hypothesis is invertible. Then $\mathbf{x} = W(\lambda) \mathbf{1}$, by $(*)$, and $\mathbf{1}^T W(\lambda) \mathbf{1} = 0$, by $(**)$. The last scalar equation, when multiplied by $\det (M^2-\lambda E)$ becomes a polynomial in $\lambda$ of degree $n$ and in general gives $n$ solutions. The constraint equation becomes $\mathbf{1}^T W(\lambda)^2 \mathbf{1} = u$. In general, two independent equations for $\lambda$ will give you no solutions.
It is also possible to examine the case when $\lambda$ is an eigenvalue of $M^2$, with eigenvector $\mathbf{e}$. Then, $(*)$ has solutions only if $\mathbf{e}^T \mathbf{1} = 0$. Lets assume that is the case and denote by $W^+(\lambda)$ the pseudo-inverse of $(M^2-\lambda E)$. Then, the general solution of $(*)$ is $\mathbf{x} = W^+(\lambda) \mathbf{1} + \mu \mathbf{e}$, which converts $(**)$ and the constraint to the system
$$ \begin{cases} \mathbf{1}^T W^+(\lambda) \mathbf{1} = 0 \\ \mathbf{1}^T W^+(\lambda)^2 \mathbf{1} + \mu^2 \mathbf{e}^T \mathbf{e} = u \end{cases} . $$
Now you have two equations for two unknowns $(\lambda,\mu)$, so in general you can expect to find non-trivial solutions.
A: 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.
