@acer 
 

1. About unevaluated `fsolve` results
   Yes, you are correct that `fsolve` may often return unevaluated if no solution exists in the specified ranges. In my case, I will test the result explicitly to ensure it is of the form `{f=..., x=...}` with numeric values. If `fsolve` returns unevaluated, I will discard that trial point. Using `eval(result,1)` is indeed a good way to avoid unnecessary evaluation at the top level.

2. Ranges of the parameters

    For my first step, I fix $q$ and $M$ at specific values.
    Then,  search for valid double layer solutions in a range of $f \in (0,1)$ (since it represents a fractional density) and $\tau > 0$ (typically $\tau \in [0.01,0.5]$).
    After that, I allow $q$ to vary slightly around its nominal value (e.g. $\pm 0.1$) to test sensitivity.
    Later, I will also scan over $M$ in a physically reasonable range, usually close to 1 (say $M \in [0.8,2]$), since the Mach number must be slightly above unity for ion-acoustic modes.

3. Acceptable range for $x$ (potential $\phi$)

    By physical constraints, eq1 contains the term $\sqrt{1-2x/M^2}$. To avoid complex values, I must require

     $$
     x \leq \frac{M^2}{2}.
     $$
    In practice, I will restrict $x$ to a small positive range, e.g. $x \in [0.001, M^2/2 - \epsilon]$.

4. Handling possible complex values**

    As you mentioned, if the square root argument or the power terms become negative, the potential turns complex. I will explicitly restrict the search ranges so that no complex evaluation occurs.

5. Existence of true roots vs. residual tolerance

    Indeed, there may be cases where the two conditions $V(x_m)=0$ and $V'(x_m)=0$ are not satisfied exactly. In those cases, I will allow a small residual tolerance, for example:

     $$
     |V(x_m)| < 10^{-6}, \quad |V'(x_m)| < 10^{-6}.
     $$
    If exact roots are not found but residuals are below tolerance, I will accept the solution as a valid double-layer candidate.

6. Possibility of optimization instead of root finding

   That is a good point. If no root is found within the ranges, I may reformulate the problem as a constrained optimization, i.e., minimize $|V(x)|+|V'(x)|$ over the admissible ranges of $f, \tau, q, M$. This could capture “near solutions” that `fsolve` misses.

@sand15 
I need to detemine those values of (qc,Tch) that satisfies f3=0. (I know that both qc and Tch are positive)
How can I do that?

Thank you for your useful hints. In continution, I have another question.
 I want to rewrite w as w2 but I can't determine the value of M0 and M1 with Maple. What should I do?(see the attached file)
s2.mw

@dharr 
Thanks for your cmment. I have two questions:

1- Why the plot has been contain 2 layers(see attached photo)?  Is it possible to remove the lower layer?
2- How can I extract data of ploted figure (M,kc and deltah) as an excel file with 3 columns?

@panke 
Thanks for your comment. V(phi) has been defined in the attached file:
Waiting for your reply.
V(phi).mw

@mmcdara 
Similar to initial W, the new W should satifies the mention conditions at phi=0 and phi_0 (e.g., W=dW/dphi=0 and d^2W/dphi^2 <0 and ....).
Also:

• Do you still want three extrema and a negative function (meaning the same pattern as your initil W)?
• Yes
• Do you already have a solution?
• Unfortunately no! Looking to find the first data.
  But I know that M should be a positive number and 0<d1<1.
  Here, the work is more complicated that the previous W.
  In fact, I am looking for a code that can change the value of each parameter)(k1, k2, d3, ...) in a given interval and see if there is a (M,d1)  that with the given parameters makes the W have three extreme like the initial W.
  ​Maybe our answer to the above question be No! there is not any values of (M,d1) for whcih we have 3 extreme points in W.
  
  The range of changes of some of the above parameters is as follows:
• k1=2..5
  k2=3..8
• s1=100-150
• s2=1..5
• s3=0.01..0.3
• d3=0.001-0.004

@mmcdara 
Thank you for your post.
To better understand your post, I tried to use your code and your comprehensive comment for a new 'W' function, but I couldn't. May you help me?(please see the attached file)
New.mw

@mmcdara 
In answer to following question:

Would you be interested in finding  the way to get, not a but several, triples T°= (phi_0°, f°, beta°) for which

• W(0, T°)=0, W'(0, T°) = 0
• W(phi_0=phi_0°, T°) = 0
• W(phi, T°) < 0, for each phi in (phi_0°, 0)
• W has only 3 extrema in the range (phi_0°, 0)

Yes. I'm interested in this. Meanwhile, depending on values of beta and f it  is possible that these 3 extreme points happen for phi>0(the right side of cordinate origin) and/or phi<0 (left side of origin) .
Also, if I consider the first extreme at phi=0 (at origin of coordinate system), your are right and we'll have 4 extreme points.
Moreover, Phi_0 and those value(s) of beta and f that satisfies the above conditions are unknown.
I hope you can help me, as you said at the end of your above comments.

@mmcdara 
Thanks for your comments. Considering a function similar to W,  what I need is to know for which value of f and beta, the function  W has a behavior similar to the plot in my uploaded file (that is, 3 consecutive maxima under the horizonal axis while w and dw/dphi=0 at phi=0) and these values can be extracted in the form of a T table.

@Carl Love 

In fact, it means that there should exist 3 extreme point beteen phi=0 and phi=phi0(similar figure in  the attached file)

@mmcdara 
Thanks.
I want to determine those vales of (f,beta) for which there exist 3 extreme point in plot of W vs phi, but we have found the location of the extreme points!

Thank you for your help.

@Rouben Rostamian  
I appreciate you for your help.
please see my new uploaded file. 
If we know [kc>3/2, kh>3/2, (M-u0b+sqrt(3)*sqrt(mu*sigma[b]))>0, M-u0b-sqrt(3)*sqrt(mu*sigma[b])>0, M-sqrt(3)*sqrt(sigma[i])>0, M+sqrt(3)*sqrt(sigma[i])>0], then for phi=0 we analytically have f1=0:

f1:=-(1-delta[h]-delta[b])-delta[h]-(delta[b]/18)*(3*mu*sqrt(3)*(2*sqrt(3*mu*sigma[b]))/(mu*sqrt(mu*sigma[b])))+(3*sqrt(3)/18)*(2*sqrt(3*sigma[i])/sqrt(sigma[i]))  # f1=diff(w,phi) at phi=0
 but with Maple f1(phi=0)<>0. 
 What is the reason for this difference? How can I get this result with Maple?
N-eval.mw

@Rouben Rostamian  
Thank you for your reply. 
The worksheet has been changed.
The uploaded worksheet was not complete, so I sent it again. Please see the new file
New-eval.mw

@Carl Love 
Thank you for your reply. 
First of all, I must say that your guess is correct.
As I wrote to you before, my function (q3) is a multi-variable function. Some variables have a fixed value (for example, u0b or mu , ...), but others, such as M or delta_h, have a range of changes and "all values" within the range must be considered, not just the first and last points of the range. My goal is to find the 'phi' or (phi)s in which the function q3 satisfies the following three conditions:
1)q3=0, dq3/dphi=0 at phi=phi_0 where 0<phi<phi_0,
and
2) q3<0 for 0<phi<phi_0.
3) q3=0, dq3/dphi=0 and d^2(q3)/d(phi)^2<0 at phi=0. 

I'm waiting for your reply.
q3.mw