I want to return {x = 0., y = 1.158748796 + 0. I} as {x = 0., y = 1.158748796 }. The solution is coming from:

soln3:= fsolve({b1, b2}, {x = 0 .. infinity, y = 0 .. infinity});

and the second solution is coming from:

soln4:= fsolve({b1, b2}, {x = -infinity ..0, y = -infinity .. 0});

See my code below

restart:

Procedure

doCalc:= proc( xi )

# Import Packages

uses ArrayTools, Student:-Calculus1, LinearAlgebra,

ListTools, RootFinding, plots, ListTools:

local gamma__1:= .1093,

alpha__3:= -0.1104e-2,

k__1:= 6*10^(-12),

d:= 0.2e-3,

theta0:= 0.0001,

eta__1:= 0.240e-1,

alpha:= 1-alpha__3^2/(gamma__1*eta__1),

c:= alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1)),

theta_init:= theta0*sin(Pi*z/d),

n:= 30,

g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,

qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,

soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,

AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar,

soln3, soln4, Imagroot1, Imagroot2;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

g:= q-(1-alpha)*tan(q)+ c*tan(q):

f:= subs(q = x+I*y, g):

b1:= evalc(Re(f)) = 0:

b2:= y-(1-alpha)*tanh(y) -(alpha__3*xi*alpha/(eta__1*(4*k__1*y^2/d^2+alpha__3*xi/eta__1)))*tanh(y) = 0:

qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):

OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);

# Compute Odd asymptote

ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));

qstarTemporary:= min(ModifiedOddAsym);

indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);

qstar2:= OddAsymptotes(indexOfqstar2);

# Compute complex roots

AreThereComplexRoots:= type(true, 'truefalse');

try

soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});

soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});

qcomplex1:= subs(soln1, x+I*y);

qcomplex2:= subs(soln2, x+I*y);

catch:

AreThereComplexRoots:= type(FAIL, 'truefalse');

end try;

# Compute the rest of the Roots

soln3:= fsolve({b1, b2}, {x = 0 .. infinity, y = 0 .. 10});

soln4:= fsolve({b1, b2}, {x = -infinity ..0, y = -infinity .. 0});

Imagroot1:=subs(soln3, I*y);

Imagroot2:= subs(soln4, I*y);

OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);

AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));

if AreThereComplexRoots

then gg:= [qcomplex1, qcomplex2,op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),

seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];

elif not AreThereComplexRoots

then gg:= [op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];

end if:

# Remove the repeated roots if any & Redefine n

qq:= MakeUnique(gg):

m:= numelems(qq):

## Return all the plots

return qq, Imagroot1,Imagroot2, p, soln3, soln4;

end proc:

ans:=[doCalc(0.06)]:

ans[5];

{x = 0., y = 1.158748796 + 0. I}

ans[6];

{x = 0., y = -1.158748796}