@acer Yes, of course, I meant that it was nearly impossible via the characteristic polynomial via the route in the OP's worksheet.

For this purpose, it would be nice to compute the dominant eigenvalues or singular values without needing to spend the effort to compute them all. Is there a way to facilitate that? 

@asmakhan What variable or quantity in your problem represents temperature? In other words, what variable do you want isoclines (isotherms) for? What are the ranges of the parameters that you want to consider?

@Stretto I was simply pointing out DataFrames as an example of a Maple object for which ...[...] is an overloaded operator that has special meaning when applied to relations (such as equalities and inequalities). You can use showstat to see how that is coded. Yes, they are like tables, but they can be indexed by relations as well as the usual indices.

@Magma Yes, I tested a 128x128, and that was the largest that I tested. It ran without problems, and didn't use any significant amount of memory: 22 Meg. It completed in 98 seconds, the final matrix was 1078 columns, and the final XOR count was 2881.

I guess that you're using 32-bit Maple also? Find the places in the code where it says datatype= integer[8]. Change that 8 to 4 and see what happens. The 4 is 4 bytes per word, or 32 bits.

Here's how to generate random test cases that are repeatable:

key:= 1:
randomize(key):
#Now generate a test case:
M:= LinearAlgebra:-RandomMatrix(2^7 $ 2, generator= rand(0..1));

Now if you do the random generation again with the same key, you'll get exactly the same M. Vary the key to other positive integers and see if you can get another example that crashes. If you find one, you can tell me the key, and I  should be able to generate it. (Although I'm not sure how many versions one can span and still maintain that continuity. Yours is 8 years older than mine.)
 

@Ramanujan The purpose of my questions about N is to get diagnostic information that will help me to make a coefficient extractor. It wasn't suggesting it as a workaround for you to use in lieu of an extractor.

@mwahab Please carefully check the accuracy of some terms with squares or higher powers of the derivatives. I may need to change the frontend argument {`+`, `*`} to {`+`, `*`, `^`}.

Regarding the nonlinearity: I wasn't saying that there was any mathematical flaw if your equations had it. I was just saying that my original code would not handle it correctly. The current code should be able to handle polynomial nonlinearities, but it still cannot handle terms such as u[x]*sin(u[x]).

@Magma By all means, suggest improvements!

But I don't see the connection between that and removing unnecessary comments. When you've earned the status of a Moderator (500 reputation points), then you can remove them. At that point, you may realize that it's unwise to do so.

@Magma Please put the original pseudocode back in the Question! Even though it had some small errors, it was far easier to understand than the C++ code. It's far better to note the mistakes in information rather than removing it entirely. I always tell my math students, "Cross out; don't erase!"

@Carl Love By just updating the columns of H that will be changed due to the changes in R, as I mentioned earlier, I've got the time down to 3 seconds for a 64x64 matrix.

Paar:= proc(M::Matrix)
option `Author: Carl Love <carl.j.love@gmail.com> 25-Oct-2019`;
local 
     R:= rtable(M, datatype= integer[8]), n:= upperbound(M)[2], xors:= add(R)-n, 
     H:= rtable(R^%T.R, datatype= integer[8]), 
     hmax, maxij, mi, mj, lastcol, Update, H_as
;
     for lastcol from 1+n do
          H_as:= rtable(antisymmetric, H);
          hmax:= max(H_as);
          if hmax <= 1 then return (R, xors) fi;
          member(hmax, H_as, 'maxij'); 
          (mi, mj):= maxij;
          R(.., lastcol):= R[.., mi] *~ R[.., mj]; #Add new column to R.
          #Update old columns of R:
          R[.., [maxij]]:= R[.., [maxij]] *~  (1 -~ R[.., [-1$2]]); 
          xors:= xors - (hmax-1);
          H(.., lastcol):= <seq(0, 2..lastcol)>; #Add column to H.
          H(lastcol, ..):= `<|>`(seq(0, 1..lastcol)); #Add row to H.
          Update:= R^%T.R[..,[maxij, -1]]; #3 new columns/rows for H.
          H[.., [maxij, -1]]:= Update; #Update H's columns. 
          H[[maxij,-1], ..]:= Update^%T #Update H's rows.
     od
end proc
:

Please verify my results very carefully for some small matrices. As I said, I'm not familiar with this material. Nor have I ever performed such a delicate "matrix surgery": appending a row and a column, inplace, on a square matrix. This surgery is the last 5 lines of the loop, and I hope that that number of lines can be reduced.

Since the number of rows/columns that need to be updated is constant, 3, regardless of the input size, this represents a reduction of the asymptotic time complexity of the original published Paar's algorithm by a whole power. Thus, this might be worth publishing.

@Magma The Maple command RowDimension is poorly named. It is the number of rows, not the dimension of the rows, which would be the number of columns. In your code, you set n to the number of rows. I believe that that's incorrect, and it should be the number of columns. Please confirm.

@mwahab The error is a result of an obscure Maple version difference that I forgot about (so obscure that I've never seen it documented). Change the table line to

Cfs:= table([terms]=~ [Cfs]):

It turns out that doing the entire matrix product H:= R^%T.R is much faster than doing just the dot products needed to get the above-the-diagonal entries. So, this version does the same 64x64 matrix in 1 minute.

Paar:= proc(M::Matrix)
local 
     R:= rtable(M, datatype= integer[8]), n:= upperbound(M)[2], xors:= add(R)-n, 
     H, hmax, maxij, lastcol
;
     for lastcol from 1+n do
          H:= rtable(antisymmetric, R^%T.R, datatype= integer[8]);
          hmax:= max(H);
          if hmax <= 1 then return R, xors fi;
          member(hmax, H, 'maxij');
          R(.., lastcol):= R[.., maxij[1]] *~ R[.., maxij[2]];
          R[.., [maxij]]:= R[.., [maxij]] *~  (1 -~ R[.., [-1$2]]);
          xors:= xors - (hmax-1)
     od
end proc
:

Further optimization is likely possible. H does not need to be entirely recomputed each iteration because only 2 columns/rows have changed and 1 needs to be added. This type of inplace updating of a growing matrix is tricky, and it'll take me a few days to get to it.

@Magma I modified the code so that it returns the updated matrix followed by the xor count. For my first example with n=64, I started with a random binary 64x64 matrix. The result was 324 columns, the xor count was 784, and it took 16 minutes to compute. If I declare the matrices R and H with datatype= integer[8], I can do it in slightly less than half that time. What times were you seeing?

Paar:= proc(M::Matrix)
local R:= copy(M), H, hmax, maxij, lastcol, n:= upperbound(M)[2], xors:= add(M)-n;
     for lastcol from n do
          H:= rtable(antisymmetric, (1..lastcol)$2, (i,j)-> R[..,i]^%T.R[..,j]);
          hmax:= max(H);
          if hmax <= 1 then return R, xors fi;
          member(hmax, H, 'maxij');
          R(.., lastcol+1):= R[.., maxij[1]] *~ R[.., maxij[2]];
          R[.., [maxij]]:= R[.., [maxij]] *~  (1 -~ R[.., [-1$2]]);
          xors:= xors - (hmax-1)
     od
end proc
:
M:= LinearAlgebra:-RandomMatrix(2^6 $ 2, generator= rand(0..1)):
CodeTools:-Usage(Paar(M));

Note that the dot product of binary vectors equals the Hamming weight of their conjunction. I've been guessing that due to processor optimizations, dot products over the integers are faster than counting the Hamming weights the other way.

@mwahab My code needs some small modifications to handle product (or "multiple") terms. And here it is:

Jxy:= [indets(eqn, specindex(op~(0,Jets)))[]]; #just an example

Eqn:= collect((lhs-rhs)(eqn), Jxy, 'distributed'):
Cfs:= frontend(coeffs, [Eqn], [{`+`, `*`}, {}], Jxy, 'terms'):
Cfs:= table([terms=~ Cfs]):

The first line will vary depending on the type of terms that you're looking for. I recommend that you use this new version regardless of whether there will be product terms. The frontend command is needed to prevent coeffs from looking inside the F, P, etc. terms.

If you can't bracket the N as requested above, could you simply provide the smallest N that you can find for which the results are undesirable and the largest N that you can find for which the results are good?