@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

@Carl Love 
Hi Carl Love,
Can you give me some advice to solve this problem that I raised with you yesterday?

@Carl Love 
Thank you very much for your comment. That error was removed, but unfortunately, my worksheet does not work well. 
Can you help me?
My problem is:
Assuming having the q3 function and the specified parameters in the attached file(kc, kh, ... delta[b]), I am looking to find a point or points (for example, phi=phi_0) where the following conditions are met for specific values of delta[h] and M in mentioned intervals and then draw q3 versus phi:
1)q3=0, dq3/dphi=0 at phi=phi_0 where 0<phi<phi_0,
and
2) q3<0 for 0<phi<phi_0.

Of course, at phi=0, the conditions q3=0, dq3/dphi=0 and d^2(q3)/d(phi)^2<0 must be met, which are met.

Sincerely.

@Rouben Rostamian  
As you mentioned, I tried to solve the above Eq. numerically, but I see two curves. why? and What should be done to draw only the linear curve?
Eq.mw

@sand15 
Thank you for your help. As I write in attached file, I want to calculate f1^4 integration, not f1. Is it possible with Maple analytically or, similar to f2, it need numerical integration?

@tomleslie 
Dear Sir,

How can I turn the three parametric equations (Eq22, Eq33, Eq44 in the attached file) into an equation in term of phi[d0] and delta[d]?
I want to get an equation similar to C1*phi[d0]+C2*delta[d]=0 where C1 and C2 are known coefficients.
Assume that u[i10], u[i20], phi[d0] and delta[d] are the only unknowns in mentioned equations and the rest of parameters (for example, alpha, delta[i],...)are known.

Equation.mw

@acer 
 

Hi,

What can be done when a question is not answered but sent back to other people??

@raj2018
Dear  tomleslie ,
Are the above parameters sufficient to solve the set of NL equations?

@tomleslie 

Sorry! 
I have been inaccurate when writing. In fact, using equations 1-4 and assuming that  alpha=0.01, delta[d]=1e-4 and the range of changes in delta[i] = 10e-5..10, I want to plot the variation of u[i01],u[i02] and phi[d0] versus delta[i].

A paragraph of the article is as follows:

"The non-linear equations (1)-(4) constitute a complete set of equations enabling the computation
of delta[e], phi[d0],u[i10], and u[i20] for different  parameters such as alpha=0.01 and delta[d]=1e-4."

@tomleslie 
Hi,
As you said, Equation 1 had a problem, and now I have solved the problem, and the number of equations is consistent with the unknowns. Please see the attachment and help solve the equations.

Regards,

H.mw

@tomleslie 

Similar to your worksheet, I tried to write my equations in the form    ” f(u[i10], u[i10], phi[d0])=delta[d]”

Equations 1 and 2 were written in this form easily, but I couldn’t do that for equation 3  due to following error:

`[Length of output exceeds limit of 1000000]`

Solve.mw

@tomleslie 
Thank you for your reply. I'll try to solve my Eqs. with your comperhensive response. 

Can I have your email address?

r.choupta@gmail.com

@tomleslie 

 

@tomleslie 
I tried to sent the article file to you, but I couldn't. However, what is your opinion about this assumption to plot the figures?

Assuming that delta[d] varies from 10^(-6) to 10^(-3), we have 3 equations and 3 unknowns u[1i0], u[2i0] and phi[d0]. Therefore, I solve the set of equations several times by a loop to find "u[i10], u[i20] and phi[d0]" for each step of the loop.

For example,

#%%%%%%%%%%%%%%%%%%%%
for delta[d] from 10^(-6) by 10^(-6) to 10^(-3) do
sys:=eval([Eq2a,Eq3,Eq4],[delta[i]=0.5,alpha=0])
fsolve(sys)
print delta[d], u[i10],u[i20],phi[d0]
end do
#%%%%%%%%%%%%%%%%%%%%


How do I do this by Maple and how can I extract data for ploting u[1i0] vs delta[d],....?
Now,

@tomleslie 
I sent you the article. Please see the attached file.(see the last paragraph in Sec. II and the curves in Sec. III)

[copied without permission, deleted by moderator]

According to the above message, How to send the article file to you?
Please, see the following link:

https://doi.org/10.1063/1.4967763

@tomleslie 

How do you think the author of the article was able to draw the relevant figure? The article only mentions these four equations and parameters that I wrote. The exact phrase of the article is:

The non-linear equations (1), (2), (3), and (4) constitute
a complete set of equations enabling the computation
of  delta[e], phi[d0], u[i10], and u[i20] vs delta[d] for different physical parameters !!