Explicit solution to a Rayleigh quotient equation

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.

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.