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

@Alfred_F A bit more, with a nice example at the end.

Download floats.mw

@Ronan I usually copy the URL here so you can give it as additional information. Then go to the "more" tab on mapleprimes and choose submit software change request. Maple employees may notice stuff here, but I would submit the SQR anyway, just in case.

I get the same one-page result as you for both 2024 and 2025; both through exporting as .pdf or printing to the Adobe print driver. (Mac .pdf generation is different I believe.)

I agree with @nm that it's a bug. The following works and is shorter:

B:=<seq(A)>;

@JaneCherrytree Thanks. My point was to simplify the expression as much as possible before integrating, and in particular getting rid of square roots, but if it is correct then we have to live with the square root. What can we assume about the signs of the variables? Are they all positive?

You have

t1 := r*cos(phi);
t22 := (1 - s)*sqrt(1 - t1*t1) + s*t2;

So  sqrt(1 - t1*t1) = sqrt(1 -r^2*cos(phi)^2) seems wrong from a dimensional point of view since r^2*cos(phi)^2 has dimensions if r^2 and is subtracted from dimensionless 1. Is it supposed to be sqrt(r^2 -r^2*cos(phi)^2), which would simplify to r*sin(phi)?