I agree that the documentation suggests that Q and rp are the same size "r". However, a look at the code using showstat(RowEchelonTransform) suggests that the length of Q (=P in the code) need not be the same as rp. P is initialized as a 1..n row vector (n=number of rows presumably) (line 7), but on line 9 only the entries from 1 to min(r,n-1) are output. Presumably r is the rank and in your case min(16,15) = 15.

Since the length of Q is the number of swaps needed to construct P, perhaps you don't need to do the last one if it is full rank, just do all the swaps in Q. I suspect this will work, but if not you could just do RowReduce as a workaround if you only want the row reduced form.

Edit: I confirmed this works for the reduced form (redflag = true). The redflag = false case gives different answers for RowEchelonTransform and RowReduce, but the nonreduced form is not unique. So I think the swap procedure is correct.

I submitted an SCR for the help page.

Download Modular.mw

Your side note question is the easiest - the spaces are missing commas. The code is something like this

P:=proc();
userinfo(1,P,`entered with args:`,args);
7;
end proc:
infolevel[P]:=1;
P(as,3,5);

which gives 

P: entered with args: as 3 5

For your main question, I don't think there is a simple answer - many internal routines are called and internal tables probably retain some information. For your first example, the message 

Solving each unknown as a function of the next ones using the order: [y(t), x(t)]

comes from `PDEtools/sdsolve`, but I don't see where the next message comes from. You could set printlevel to a high value to track it down, but even if you did, the solution (perhaps forget on a deep routine?) would have to be different for different odes.

(I encountered the same problem with SolveTools:-Polynomials system and infolevel[solve].)

output=basis is for exact solutions; earlier you had integer parameters and an analytical solution in terms of Heun functions was possible (not sure if this is the same equation). So you could try leaving those parameters as variables, find the exact solution, and then plug the values of the parameters into the solutions. I gave up on this after some time computing, but it might ultimately work. (The general solution is then a*(solution 1)+b*(solution 2), where a and b depend on the initial conditions.) Otherwise, for a numerical solution you need to choose initial conditions, as you succeeded to do in an earlier case.

My interpretation is that convert(...,unit_free) only removes the units from a product of a number and unit. So it removed the 1/um^2 but not the unit in the argument of the exponential ("data structure with units embedded in it").

subsindets(expr, anything, convert, unit_free)

will remove all the units.

You had diff(E,psi), which evaluates to 0 instead of diff(E(psi),psi). Note also gamma has a special meaning unless you use "local gamma". Changes in red. I didn't proceed beyond the mess you now get for fin1.

Note: Maple's JacobiSN is spelled with J not j, and has two arguments, so you will have to fix that.

0123.mw

@MaPal93 sent me some notes about how argmin was derived, so I could proceed with a matrix-based method, similar to in the earlier thread. This can probably be refined for the second part of the problem, but I used it in the maximization part to avoid ever creating a problematic denominator, thus avoiding the issue. For the calibration problem, there are a few restrictions on lambda__2, then solve a quadratic for lambda__3, then a linear equation for lambda__1. So finding conditions for positive solutions should not be intractable, though perhaps tedious.

MatrixNew2.mw

Edit: There are no positive solutions for the calibration problem (meaning the conjecture is not correct), and a more detailed analysis is given below.

Unless there is some dramatic cancelling of terms, then probably not. A 32x32 determinant has 32! terms. Now 32!/16! = 12576278705767096320000 = 1.26 x 10^22, which is how many times longer the calculation for a 32x32 determinant will be over a 16x16 one as an estimate.

On the other hand, your matrix has a lot of structure, which you presumably know about in detail. Did the 16x16 one have 16! terms or was it much simpler? Do you have a guess as to what the answer is, or what form it takes? Perhaps it can be divided into simpler tasks (Determinants of submatrices, perhaps). Do you really want this determinant, or is it just a step to something else that can be calculated without it (solutions to a linear system for example)? Do you just need some information about it? - zero or non-zero perhaps? Do you know anything about the signs of the variables?

Without some guidance about the nature of the problem, which might suggest a clever solution, I'd say you are out of luck. 

This partly duplicates @mmcdara's analysis (I voted up!) that used numerical values for the parameters, but with the analytical solutions, and shows how to interpret the PolynomialSystem answer and do the backsubstitution (part of your second question). In this case solve is simpler (and presumably the same as if backsubstitute=true was chosen, but I didn't check), but had there been restrictions on lambda__2 and lambda__3, then fishing out the required solutions from all of them probably would have been easier starting from the PolynomialSystem, backsubstitute=false form, as was the case in your earlier problem.

Download solvers.mw

Looks like a bug. As you probably know, you can set the display precision for the whole worksheet from the menu Tools -> Options -> Pecision.

The divide symbol is not on my keyboad. Yours appears as "÷" in the string. Assuming yours is always encoded in the same way, then it can be substituted with "/" using StringTools:-Subs. I needed to use 1-D entry here.

PEMDASTest.mw

It's not obvious which variables to eliminate, but this works. 

Download elim.mw

For a numerical solution the limit boundary condition can't be used. If you want to approximate infinity by by a large number you can use, say, U[2,n](20)=0. But then you have boundary conditions at -Pi, 0 and 20; the solver needs just two boundary locations, so I replaced it with a boundary condition at 0, which you will need to modify to what you want. You need also to replace x[01] with X[01] to avoid confusion with the simple variable x. Then it is possible to get a solution.

Download Thesis_(1).mw

So here is a workaround,  which just loops as @sursumCorda suggested.

Download PrincipalPart.mw

Not sure how you want to label them, but here are a couple of simple possibilities, using the caption as a label.

Download Plotlabels.mw

@lcz For IsSubgraphIsomorphic, Maple uses a constraint algoritham and SAT solver, perhaps this algorithm? doi 10.1007/s10601-009-9074-3 or here? doi: 10.1016/j.artint.2010.05.002, but I'm not sure there is an equivalent for the induced case.

The VF2 algorithm and its (much) improved variants VF2plus, VF2++ and VF3, seem to be widely used. Here I only tried VF2, and didn't try to optimize the data structures (mainly sets straight from the paper); probably moving to one of the improved algorithms would be the next step. I didn't check it on a large number of cases; perhaps you have some other cases to test it.

[Edit2: Updated version here now implements some parts of the VF2++ algorithm, and removes redundancy in search tree]

It seems fast to find matches if there are some, but of course is slower to show there are no matches. Hope this is useful.

[Edit - full VF2++ below]

Download VF2conndegAlgorithm4.mw