@sand15 

Thank you for your answer. As you showed, there may be several values of $f$, $\tau$, and $\phi$ that satisfy both Eq. (1) and Eq. (2). Among these possible values, I am looking for those that also satisfy the following conditions:

$$
\text{Eq1}(x=0) = 0, \quad 
\text{Eq1}(x=x_1) = 0, \quad 
\text{Eq2}(x=0) = 0, \quad 
\text{Eq2}(x=x_1) = 0,
$$

and

$$
\frac{d^2}{dx^2} \text{Eq1} < 0 \quad \text{at } x=0 \text{ and } x=x_1,
$$

where $x=0$ and $x=x_1$ are the two extreme points of Eq. (1).

I know that for $q=0.9$, $M=0.9973$, $f=0.01$, and $\tau=0.09$, all the above conditions are satisfied, but I do not know how to implement this procedure in Maple.

@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