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

@mwahab What you want just involves slight modifications of the indets command. I think that your question #1 means that you want the 4th derivatives of u and v. Then

Jets:= {u,v}(x,y,z,t);
Jxy:= indets(eqn, myindex~(op~(0, Jets), identical(op~(Jets)[])$4));

The only difference is $4 instead of $2. If you want all derivatives of u and v, then

Jxy:= indets(eqn, specindex(op~(0,Jets)));

@mwahab What do you mean by the phrase "and their multiples", which you've used several times? Please show a specific example where a derivative has "multiples".

@Ramanujan Could you do some manual binary search experimentation to find the largest N for which you get good results, or conversely, the smallest N for which you get bad results?

Your series command is meaningless unless eta is expressed in terms of q. In order to help you, we need to see that expression.

@Magma Change the + in the exponent to %T

@Carl Love Here's a better way to automatically generate the numeric evaluation procedure for a first-order recurrence. This doesn't need a remember table because it always makes only one reference to the previous member of the sequence.

Makeproc:= proc(Reqn::`=`, IC::function(integer)=complexcons, Depvar::function(name))
local 
     A:= op(0, Depvar), n:= op(Depvar), An,
     P:= subs(
          _R= unapply(subs(A(n-1)= An, solve(Reqn, A(n))), An),
          (n::integer)-> _R(thisproc(n-1))
     )
; 
     assign(subs(A= P, IC));
     eval(P)
end proc
: 
SP:= Makeproc(S(n)-S(n-1) = 2*S(n)^2/(2*S(n)-1), S(1)=1, S(n)):
SP(100) - SP(99); 
SP(10000) - SP(9999);

@Magma Here is code retrofitted to not use max[index], so I believe it'll work in Maple 15:

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

@Magma This code has been retrofitted for Maple 2015. I don't actually have Maple 2015 to test that, so please test it.

Paar:= proc(M::Matrix)
local H, R:= copy(M), maxij, newcol;
     do
          H:= rtable(antisymmetric, R^+.R);
          maxij:= [max[index](H)];
          if H[maxij[]] <= 1 then return R fi;
          newcol:= R[.., maxij[1]] *~ R[.., maxij[2]];
          R:= <R | newcol>;
          R[.., maxij]:= R[.., maxij] *~  (1 -~ <newcol|newcol>)
     od
end proc:

@max125 I don't know why multiplying by the denominator makes rsolve able to solve it. It seems like a very old and not-well-maintained part of Maple, if it's maintained at all.

Your first solution attempt failed because in order for an rsolve command to make any sense, n must be a variable, not a number.

Your second attempt failed because the procedure returned by rsolve(..., makeproc) is a simple recursive procedure that can only work when n is a number, not a variable. The procedure returned by makeproc (in this, nonlinear, case) is essentially equivalent to the one returned by the following technique, which can be applied to any first-order recurrence:

restart:
Reqn:= S(n)-S(n-1) = 2*S(n)^2/(2*S(n)-1):
IC:= S(1)=1:
SP:= subs(
     _R= subs(S= thisproc, solve(Reqn, S(n))), 
     proc(n::posint) option remember; _R end proc
);
assign(subs(S= SP, IC)); 
SP := proc(n::posint)
option remember; 
   thisproc(n-1)/(1+2*thisproc(n-1)) 
end proc
SP(100) - SP(99);
                              -2  
                             -----
                             39203

Experts' note: Note that it's possible to use thisproc as a name outside a procedure and have it assume its usual role when subs'd into a procedure. This was a refreshing surprise.

There are two problems with your posted algorithm: The first is that mxcoli and mxcolj should be maxcoli and maxcolj. The second is, I believe, that !(newcol & maxcoli) should be (!newcol) & maxcoli, and likewise for j. I looked up the algorithm online, and I found a source to confirm this; but you should doublecheck because I'm not familiar with this material.

@Carl Love I prefer Acer's use of the undocumented command PiecewiseTools:-Support over my use of "magic", although documented, numbers.

@mwahab If I load your worksheet (making no modifications whatsover) and press !!! from the toolbar (Execute Entire Worksheet), then I get the expected results.

It's a good idea to start any worksheet with a restart command, alone in its own execution group.