@mary120 I cannot read the equations and conditions properly. But you said "For some values of parameteres, the output curves were very similar to what I expected". I took that to mean you have some suitable ouput. I would like to see that, because it would clarify the questions I asked and your reply does not help.

Did you try the code provided here? (Hint: you may have to simplify the result.)

@mary120 @mmcdara has suggested you were really wanting a different differential equation than the one you specified in your worksheet, and since you have changed his answer to best I'm assuming that is correct? But now you are saying you can get the correct answers with my worksheet, which uses one branch of the original ode. So is that the right ode and the right branch? Which was the correct initial condition? I suggest you upload a corrected worksheet in which your parameters that do give close to the required solution are given, so we can see the answers to these questions.

The answer to your immediate question is that to get V(phi) values down to phi=0, you need to solve the de all the way out to large enough x values that phi is close to zero, i.e., choose a wider range in

From your graph of phi vs x at initial condition phi(0) = 0.022, it looks as though going to x = 250 is about right.

@mmcdara You are pobably right, but the OP titled the question with "first oder ode", which I took at face value.

@salim-barzani OK, I understand. Here's how to automate the substitutions. But it is not a solution. Each of your pdes needs some individual attention, to find independent equations. For most cases, the expx in f are the same exps that survive after substituting into the pde, but if thay get messed up, then one needs to work harder, as here. I don't know the general solution to this issue. It may be that simplification with side relations can work; I'll give this some more thought.

Dr.D-pde.mw

I guess my alien powers have deserted me!

@salim-barzani If you only want the variables, and not the more complicated expression, use indets(something, name).

I think this is what you are asking, but if it is not, please be more specific about what you want.

@mmcdara I'm not sure of the point you are making. If it is just that you have to be careful of branch cut issues, I agree. In this case the solve strategy I used to answer the OP works - if both solutions had been real there would have been less of a problem.

(Strangely, the worksheet is not displaying, in a way I haven't seen before.)

Download DE_solve.mw

Perhaps the simplified solution is not correct? When I substitute in, I don't find everything zero. The nonzero equation is very similar to the rifsimp expression, but with 4 instead of 8.

Download demo_question.mw

@mmcdara But you have ode := diff(phi(x), x) : sqrt(1-phi(x)): instead of ode := diff(phi(x), x) = sqrt(1-phi(x)): Then the numeric one does agree with one of the two solutions.

@mmcdara Very nice, vote up. You probably meant B[j[]] not B[j] (as for theta).

@salim-barzani If you need to, you can use assuming n::posint to assume n is a positive integer, assuming n::positive or assuming n>0 for n positive.

Does simplify give you what you want?

In general, Maple will not assume n is a positive integer unless you tell it. (Not needed here, though.)

Internally z is Float(16,-3); op(0,z) gives Float and op(z) gives 16,-3. lprint(z) gives .16e-1, which I assume is some sort of canonical form.

So the 0.016 is just for display. While it might be possible to recreate that display routine, probably sprinf gives you the ultimate in control for this.

@scottwr For inner nested loops, use print for the things you want to see. So in this case just use print(x*y) instead of x*y. 

@mmcdara That was my comment. In the paper they condense products of subscripted symbols into single subscripted symbols, where the new subscripts are the set union of the subscripts of the factors. After looking further at the paper, I think it is a bad idea to do this since one loses information.