@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.

@Carl Love True, but the standard way to abbreviate input for anything nontrivial would be with diff_table. For example

declare(u(x));
U:=diff_table(u(x));

then U[x,x] is used to abbreviate diff(u(x),x,x) and U[] for u(x).

@FZ Sorry, I don't understand what they are doing.

@FZ Please expand the image. As far as I can see, they use the first integral, which didn't work for your case unless you know the form for the integrating factor.

@FZ Your equation has a derivative on the rhs, unlike your textbook example, which can be done like this:

Download Hamiltonian3.mw