@acer 

Oh I see

@acer Okay thank you. I quite understand my project. I tried to do the same work in Matlab. I was having a different result which made me check my equations and I discovered there is a change of sign in the g equation. I used maple to check the analytical solutions and Matlab for numerical solutions. That was when I discovered the error. The only problem I have is the fact that when I change the sign from (-) to (+) in g equation maple only returns the real roots instead of the real and imaginary parts. When I solve it with Matlab, I discovered the function is purely imaginary and maple is not finding those roots. 

@acer 

While I'm passing judgment because the function r(q) is a complex function that is purely imaginary and Maple can't find all the roots. If there is any way I can buyout the process, I will appreciate it. 

Secondly, why is maple taking a long time to compute this: (gridji,rngj,rngi):=[81,81],0..3,-2.0..-0.0;

It took me more than 24 hours to do this (81 by 81 grids instead of 5 by 5) on my machine. Is there any way to shorten the time?

@acer 

Thank you very much. I didn't bother to send the version I corrected the x-axis because I will also need the negative region.. 

@acer 

I'm seriously having an issue with maple 18. I used this: r:=LinearSolve(A,B). But it returns undefined values if used that, get breaks in my contour plot. see the attached plot, which differs from the one you have.

@acer 

Thank you for this:  r:=LinearSolve(A,B).

Thank you I have been able to sort it. 

@tomleslie 

Thank you

@Carl Love 

Okay, thank you.

@tomleslie 

Sorry this  is the right code

 

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, w, 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

                       qq := [2.106333379+.6286420119*I, 2.106333379-.6286420119*I, 4.654463885, 7.843624703, 10.99193295,14.13546782, 17.27782732, 20.41978346, 23.56157073, 26.70327712, 29.84494078, 32.98658013,         36.12820481, 39.26982019, 42.41142944, 45.55303453, 48.69463669, 51.83623675, 54.97783528,                     58.11943264, 61.26102914, 64.40262495, 67.54422024, 70.68581510, 73.82740963,                                 76.96900389, 80.11059792, 83.25219177, 86.39378546, 89.53537903, 92.67697249];
                        m:= numelems(qq):

                        for i from 1 to m do 
                        w[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2-alpha__3*xi/eta__1);
                        end do;

## Return all the plots
  return w;
  end proc:

ans:=[doCalc(0.06)]:
ans[1];

@tomleslie 

Thanks Tomleslie

@Carl Love 

Thank you so much

@Kitonum 

Thank you for this. But how do delete the zero root in the list of the roots? Or is there any way I can the roots that are not zeros

@Carl Love 

Thank you. Please, how do I delete zero entries from the list?  I don't want to use this method. Let's 

assume I don't the index of the zero. How do I delete it automatically without specifying the index. 

  soln=   [-0.2381948477, -0.0001434221912, 0., 0.0001434221912, 0.2381948477]

subsop(3 = NULL, soln);

@tomleslie

Hi Tomleslie,  Thank you so much for this. 

This has actually given me a clue of what is happening.  But my major problem is that it seems I have some missing roots of the function g   in the attached maple file. From the plot below it seems we have a missing root between  (1.570796327, 2.144761060)  Asymptotes. This has been a very serious challenge for me. How can I find the missing?

New_try_Ijuptil.mw

@dharr 

Thank you very much. I can pick it from here.