Specify the range of solution while using NSolve

Question

Ν = 6; M = Ν/2; prec = 20; c = 0.3 + 0.002*I; workprec = 200; 
solu = Flatten[(Λ[#1] & ) /@ Range[M] /. 
     NSolve[Table[-Product[(Λ[δ] - Λ[γ] + I*c)/(Λ[δ] - Λ[γ] - I*c), {δ, 1, M}] == 
        ((Λ[γ] - 1 - I*(c/2))/(Λ[γ] - 1 + I*(c/2)))^Ν*((Λ[γ] - I*(c/2))/
          (Λ[γ] + I*(c/2))), {γ, M}], Table[Λ[γ], {γ, M}], 
      WorkingPrecision -> workprec]]; 
ListPlot[Transpose[{Re[solu], Im[solu]}]]

While plotting the solutions to these equations, one gets

You can see that there are pair solutions of the form $\Re[\Lambda]\pm \Im[\Lambda]$ i.e. the ones which have the same real parts but the imaginary part is positive and negative (nearly equal). But there also solutions that do not come in such pairs.

I am interested in only those solutions which do not come in pairs. One way I can do is find all solutions and select the ones that do not come in pair but the problem is that it is very slow to do it that way. I have to eventually run this program for $N=100$ but it is struggling at $N=6$.

My question is:

How do I make the program search the solutions which have $|\Im(\Lambda)|<0.002$. The goal is to search the solutions only in this narrow window (where most of the non-paired solution lie) so that the program is faster. I do not care about missing some solutions or having pair solutions-- all I am looking for is making the program fast by not searching too many pair solutions so that I can scale the program to $N=100$ or above.

If you run this program, it might take quite a bit of time to get result

Update 1:

Following Michael's comment, I tried

Ν = 6; M = Ν/2; prec = 20; c = 0.3 + 0.002*I; workprec = 200; 
solu = Flatten[(Λ[#1] & ) /@ Range[M] /. 
     NSolve[Table[-Product[(Λ[δ] - Λ[γ] + I*c)/(Λ[δ] - Λ[γ] - I*c), {δ, 1, M}] == 
        ((Λ[γ] - 1 - I*(c/2))/(Λ[γ] - 1 + I*(c/2)))^Ν*((Λ[γ] - I*(c/2))/
          (Λ[γ] + I*(c/2)))&& Abs[Im[Λ[γ]]]<0.002, {γ, M}], Table[Λ[γ], {γ, M}], 
      WorkingPrecision -> workprec]]; 
ListPlot[Transpose[{Re[solu], Im[solu]}]]

This indeed printed the solutions that I wanted but this is not any faster. I could not get even N=8 in a reasonable time.

Answer 1

If you're looking for real solutions only, you may want to explicitly expand your variables into their real and imaginary parts, e.g. doing

equations /. \[CapitalLambda][i_] :> \[CapitalLambda]re[i] + I \[CapitalLambda]im[i]

Then you can separate each of your equations into two, one for the real part of the equation and another for the imaginary part (to do this very explicitly, multiply and divide each term by the conjugate of its denominator to bring all the I's upstairs).

Finally, setting \[CapitalLambda]im[i_]:>0 will leave you with N equations over N real variables, i.e. the \[CapitalLambda]re's, and because I won't appear anywhere NSolve and FindRoot will automatically look for real solutions.

WARNING: Your equations seem very much to be Bethe ansatz equations, in which case the number of solutions is related to the Hilbert space dimension, which usually grows exponentially with N. If this is the case, no matter how much you massage your equations you'll never be able to efficiently solve the full system of equations (what is attempted by NSolve) for N=100, and the way to go is probably to look into Bethe strings and the thermodynamic Bethe ansatz... though you may still find some solutions to the equations using FindRoot and some cleverly chosen starting points.

EDIT: Sorry, I hadn't noticed you had a complex coupling. In any case, to reproduce Fig. 2 from the paper cited in the comment below first define the Bethe equations

k[n_Integer, m_Integer][a_Integer] := (n - m)/2 - (a - 1);
eqs[n_Integer, m_Integer][\[Rho]0J_Complex] := 
 Table[n \[Theta][\[Lambda][a]] - 
   2 \[Pi] k[n, m][a] + \[Theta][\[Lambda][a] + 
     1/-Tan[\[Pi] \[Rho]0J]] - 
   Sum[\[Theta][(\[Lambda][a] - \[Lambda][b])/2], {b, m}], {a, m}]

Then linearize the system taking 2ArcTan[2x] -> 4x, i.e.

eqs[10, 5][-.3 + .1 I] /. \[Theta] -> (4 # &) // Expand
(*{(-13.1936 + 1.77261 I) + 36 \[Lambda][1] + 2 \[Lambda][2] + 
  2 \[Lambda][3] + 2 \[Lambda][4] + 
  2 \[Lambda][5], (-6.9104 + 1.77261 I) + 2 \[Lambda][1] + 
  36 \[Lambda][2] + 2 \[Lambda][3] + 2 \[Lambda][4] + 
  2 \[Lambda][5], (-0.627215 + 1.77261 I) + 2 \[Lambda][1] + 
  2 \[Lambda][2] + 36 \[Lambda][3] + 2 \[Lambda][4] + 
  2 \[Lambda][5], (5.65597 + 1.77261 I) + 2 \[Lambda][1] + 
  2 \[Lambda][2] + 2 \[Lambda][3] + 36 \[Lambda][4] + 
  2 \[Lambda][5], (11.9392 + 1.77261 I) + 2 \[Lambda][1] + 
  2 \[Lambda][2] + 2 \[Lambda][3] + 2 \[Lambda][4] + 36 \[Lambda][5]}*)

It's not hard to see that the coefficient of \[Lambda][j] in equation i is 2 + KroneckerDelta[i, j] (3 n + 4), so that we can solve the linear system as follows

linsol[n_Integer][\[Rho]0J_Complex] := 
  LinearSolve[
   Table[2 + KroneckerDelta[i, j] (3 n + 4), {i, n/2}, {j, 
     n/2}], -eqs[n, 
        n/2][\[Rho]0J] /. \[Theta] -> (4 # &) /. \[Lambda] -> (0 &)];

We can now use this solution as a starting point for FindRoot on the full non-linear system,

sol[n_Integer][\[Rho]0J_Complex] :=
 Block[{m = n/2, init, sol},
  init = linsol[n][\[Rho]0J];
  sol = FindRoot[eqs[n, m][\[Rho]0J] /. \[Theta] -> (2 ArcTan[2 #] &),
     Transpose@{\[Lambda] /@ Range@m, init}];
  \[Lambda] /@ Range@m /. sol
  ]

Then

ComplexListPlot[sol[200][N[-(3/10) + 1/10 I, 30]], PlotRange -> {{-1.5, 1.5}, {-.05, .01}}, AspectRatio -> 1/GoldenRatio]

produces Note however that this is still one particular solution of the Bethe equations, and you could have many others by making different choices for the quantum numbers...