May I have a answer sooner?

@acer 
The assumptions are as follows:
kappa>3/2
All variables are positive and real (phi0 is the only negative variable)
deltab and deltae are les than 1.
Is it possible to simpify the mentiond relation by this assumptions?

@dharr 
Thank you for your comments. 
If I want to simplify f in the form (M^6-const)/(M^2*(M^6-const)), how can I find the "const"s  in the numerator and denominator?

@dharr 
About your question: "But you didn't explain the significance of x=0 so it is not clear how the x axis is defined" I should say that

phi=phimax/2 at x=0 and, F(phi)=0 and dF(phi)/dphi=0 at phi=phimax

(Please, see my attached 34.mw file in one of my above reply to observe the favorite phi-x curve)

Also, I tried to attached my worksheet in that I tried to plot phi-x curve in symmetric interval, but the mapleprimes.com does not support attached file since last night. So, I had to copy it here.

restart;
F1:=phi-> 3.924999-0.24999e-1/sqrt(1-2*phi)-3.900/((1-(1/6)*phi)^(3/2))-1.648094618*10^(-14)*sqrt(3)*sqrt(1836)*(((1.3972+sqrt(3)*sqrt(1836))^2+3672*phi)^(3/2)-(1.3972+sqrt(3)*sqrt(1836))^3-((1.3972-sqrt(3)*sqrt(1836))^2+3672*phi)^(3/2)+((1.3972-sqrt(3)*sqrt(1836))^2)^(3/2))-(1/18)*sqrt(3)*sqrt(300)*(((1.4472+(1/300)*sqrt(3)*sqrt(300))^2-2*phi)^(3/2)-(1.4472+(1/300)*sqrt(3)*sqrt(300))^3-((1.4472-(1/300)*sqrt(3)*sqrt(300))^2-2*phi)^(3/2)+(1.4472-(1/300)*sqrt(3)*sqrt(300))^3)*(sqrt(3)*sqrt(300)*(((1.4472+(1/300)*sqrt(3)*sqrt(300))^2-2*phi)^(3/2)-(1.4472+(1/300)*sqrt(3)*sqrt(300))^3-((1.4472-(1/300)*sqrt(3)*sqrt(300))^2-2*phi)^(3/2)+(1.4472-(1/300)*sqrt(3)*sqrt(300))^3));
#
F2:=phi->3.927999-0.23999e-1/sqrt(1-2*phi)-3.904/((1-(1/6)*phi)^(3/2))-1.648094618*10^(-14)*sqrt(3)*sqrt(1836)*(((1.401517+sqrt(3)*sqrt(1836))^2+3672*phi)^(3/2)-(1.401517+sqrt(3)*sqrt(1836))^3-((1.401517-sqrt(3)*sqrt(1836))^2+3672*phi)^(3/2)+((1.401517-sqrt(3)*sqrt(1836))^2)^(3/2))-(1/18)*sqrt(3)*sqrt(300)*(((1.451517+(1/300)*sqrt(3)*sqrt(300))^2-2*phi)^(3/2)-(1.451517+(1/300)*sqrt(3)*sqrt(300))^3-((1.451517-(1/300)*sqrt(3)*sqrt(300))^2-2*phi)^(3/2)+(1.451517-(1/300)*sqrt(3)*sqrt(300))^3)*(sqrt(3)*sqrt(300)*(((1.451517+(1/300)*sqrt(3)*sqrt(300))^2-2*phi)^(3/2)-(1.451517+(1/300)*sqrt(3)*sqrt(300))^3-((1.451517-(1/300)*sqrt(3)*sqrt(300))^2-2*phi)^(3/2)+(1.451517-(1/300)*sqrt(3)*sqrt(300))^3));

plot([F1(phi),F2(phi)],phi=0.0..0.14);
phimax1:=fsolve(D(F1)(phi)=0,phi=0.08..0.14);
phimax2:=fsolve(D(F2)(phi)=0,phi=0.08..0.14);
integrand1:=1/sqrt(-2*F1(phi)):
integrand2:=1/sqrt(-2*F2(phi)):
ode1:=diff(x(phi),phi)=-integrand1:
ode2:=diff(x(phi),phi)=-integrand2:
ans1:=dsolve({ode1,x(phimax1/2)=0},numeric):
ans2:=dsolve({ode2,x(phimax2/2)=0},numeric):
plots:-odeplot([ans1,ans2],[x(phi),phi],[1e-2..phimax1,1e-2..phimax2]);

xmax1:=eval(x(phi),ans1(phimax1));
xmax2:=eval(x(phi),ans2(phimax2));
p1:=plots:-odeplot(ans1,[x(phi)-xmax1/2,phi],-1.4*phimax1..1.4*phimax1):
p2:=plots:-odeplot(ans2,[x(phi)-xmax2/2,phi],-1.4*phimax2..1.4*phimax2):
display(p1,p2);
M1:=plottools:-getdata(p1)[3];
M2:=plottools:-getdata(p2)[3];
 

@dharr 
If I have two different F(phi) functions (for example, see F1(phi) and F2(phi) in attached file):
1- How can solve these two ODE(s) and  plot phi vs x in one phi-x coordinate?

2- How can I plot phi-x in a symmetric interval of x (for example, from x=-100..100)
3- How can I export data of plot of these ODEs in the form of two distinct ascii table? (I couldn’t export data of the last figure in your last uploaded worksheet in the form of a ascii table)

Two_function.mw

@mary120 
I'll try to do your comment

@dharr 
Thanks for your answer. It was very useful.
How can I export the data of the phi-x plots in your attached file in  the form of an ASCII file or a two-column .txt or .dat file? 

@dharr 
Considering the problem as an ODE is a good suggestion. In fact, assuming (dphi/dx)**2+2F(phi)=0, I want to plot phi versus x.
(From your first plot, it is clear that F(phi) and dF(phi)/dx =0 at both phi=0 and nearly phi=0.09)
For an assumed F(phi), I expect the phi-x curve should be as follows(see the second Fig. in the attached file), but despite having F(phi), I could not draw the curve of phi  vs. x.

34.mw

@Rouben Rostamian  
Thanks for your comment. There was a mistake in my uploaded worksheet and I correted it (Please, see the following attached file).

Integral-New.mw

@tomleslie 
Great!
Tahnks a lot!

@Rouben Rostamian  
Thank you very much

@Preben Alsholm 
Thanks for your comprehensive comments.
If:
V:=-n^2*((T[e]+T[i])/T[e])*((((n)*(1+1/(2*delta[p]))))-ln(n)-ln(n/(2*delta[p]))-1-1/(4*delta[p])-delta[p]);
Can I expect to come up with an answer similar to soliton waves?
How can I plot the curves for different values of delta[p] in the same frame?(for example, delta[p]=0.7,0.9)

 

@Preben Alsholm

Thanks.

You are right about giving up the IC. 

If I change the IC(s) to n(x)=0 when x=+/- infinity, what happen for curves corresponding to d[p] =0.2, 0.25? 

@Preben Alsholm 

Great!!
How I can delete the shift of curves toward the RHS in such a way that the maximum of all curves locate on the y axis(n(x))?

@Preben Alsholm 
Thank you.
I try to plot RES1 for two different deta[p], namely 0.2, 0.25, by -display command, but I cauldn't. How can I do this?