Numerical method to find the roots of the expected value of a spin 2j state

Currently I am working with finding the solutions for the following problem:
I have a unit sphere in which I have n points defined by their polar and azimuthal angles: $$\theta_n , \phi_n$$. I then do the inverse stereographic projection to find the corresponding points in the complex plane: $$z = \frac{x_1+ix_2}{1-x_3}$$ where $$x_1=cos(\theta)sin(\phi)$$, $$x_2=sin(\theta)sin(\phi)$$, $$x_3=cos(\phi)$$, and then use this complex numbers as roots for a polynomial of degree n.

After this I use said polynomial to construct a quantum state $$\lvert \psi\rangle=c_1\lvert+j\rangle$$+...+$$c_{2j}\lvert-j\rangle$$, where the coefficients multipliying the base states are the coefficients of the polynomial constructed before.

I am then trying to get the roots for the following function: $$\langle S_x^2\rangle+\langle S_y^2\rangle+\langle S_z^2\rangle$$, this is going to be a function of 2*n variables and I have no idea how to approach this with numerical methods.
I would be really thankful if someone with more experience in the subject could help me with just and idea of whatever I would have to research in order to approach this.

• I wonder if it helps to redo the function in terms of $S_z$ and the raising and lowering operators $S_x\pm iS_y$. – Carl Brannen Apr 8 at 7:19
• can you please clarify what you mean by $\langle S_x^2\rangle$. Is $S_x^2$ sanwitched between two states $|\psi \rangle$ given in the question?. – David Bar Moshe Apr 10 at 9:32