Interesting. Nested radicals frequently give me problems. It is interesting that solve thinks they are the same.

expr:=sqrt(x+2*sqrt(x-1));
expr2:=sqrt(x-1)+1;
solve(expr2-expr,x);

gives x.

Probably you can use polygon from the plottools package for this, but if you upload your worksheet using the green up-arrow in the Mapleprimes editor, someone will help you.

@salim-barzani You didn't attempt to use Diff - all your Q's are zero because your B1 are zero.

@salim-barzani Does Diff work for you?

Download Diff.mw

@salim-barzani I didn't understand that. Can you give a very simple example of where conjugate doesn't allow you to do something which you want to do?

@salim-barzani You say "for diff(conjugate(u(x,t)),x) i use v[i,x]", but diff(conjugate(u(x,t)),x) is a function of t but v[i,x] is not. So for example the Laplace transform wrt t of v[i,x] will not be correct. If you fix these things up, I'll take another look.

I don't exactly understand what you want, but here are some remarks:

1. The calculation needs to be done iteratively (I think the paper uses recursively where I would use iteratively). So you need to work out all the n=0 things before the n=1 and so on. You can have one big loop but you cannot work out, say, A[0]..A[4], before you work out u[0](x,t)..u[4](x,t).

2. For debugging at least I would not use "declare"; it hides what is really happening. For example u[0](x,t) will appear as u[0], but so will u[0]. When you define B[0], you use u[0] but I think this has to be entered as u[0](x,t).

3. You ask about mixing u[i] with u[i,x]. It seems you should never have u[i.x] - only ever one subscript. The display of derivatives as subscripts doesn't mean they are actually subscripts (see point 2). The paper says u(0,x) is an initial condition, but that should be u(x,0) since you have functions of (x,t). So you can have u[2](x,0) but never anything like v[0,x] (which should be diff(v[0](x,t),x) or perhaps it means an initial condition v(x,0) but for which there presumably is a missing subscript.)

4. It is bad practice to use both u(x,t) and u[0](x,t) since they are not independent of each other; instead use U(x,t) and u[0](x,t). It may work if you are not going to assign to u, but things are not what they seem - try u[0](x,t):=g; then eval(u);

@Christopher2222 Well, maybe the symbols were generated by some other calculation which might or might have made them co-linear, but I agree that getting the assumptions right makes using the geometry package frustrating. The hints about what to assume are not always as clear as in this case.

@rockyicer  You can make your permutation matrix for the permutation [1,2,3,4] -> [2,1,4,3] using.

P := IdentityMatrix(4)[[2, 1, 4, 3]]

or you can just do the permutations of the rows and columns by

Y4_R := Y4[[2, 1, 4, 3], [2, 1, 4, 3]];

@acer Nice. I thought I had tried that, but I guess not.

@emendes assuming real is implied by beta>0. You need rho>1. It works as modified below. Notice that solve may ignore the assumptions unless you add the "useassumptions" option for solve. However, the solution is in a different form suggesting that it does take some note of them. If you do use the "useassumptions" option, it fails to solve in a reasonable time.

So for this reason I only use assuming clauses after solve and never assume before.

Download solve.mw

@Saalehorizontale It's a bug. I put it in a loop and for some reason the 9th one gives a problem (after verifying that index = 2 is correct). I have submitted a bug report.

RootOf_with_Error.mw

@Saalehorizontale Sorry, I mislead you. For a RootOf without an index, the index=1 is implied for evalf or other contexts. But you are right that it can be any of the roots of the minimal polynomial, which was just a polynomial until I turned it into a RootOf. You can check through the roots of the minimal polynomial with an exact calculation without having to worry about numerical issues; see attached.

RootOf.mw

@Carl Love @acer Thanks for the clarification; I've had that misconception for a long time. But why not have the output of p:=plot(...); the same as the actual plot? In most cases p:=foo(..); p; would give duplicate outputs.

@Saalehorizontale If not specified, it is the index = 1 root. I added a line to my answer above to check that. The order is given in the ?RootOf,indexed help page