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

@mwahab The code that I posted presumes that all the same initial commands from your worksheet--the commands that came before your definition of wave--have been executed, but your code in the Reply immediately above lacks the commands

vars := x, y, u[], u[1], u[2], u[3], v[], v[1], v[2], v[3];
PDEtools[declare]((F, P, Q)(vars));

If I include these two commands, then Cfs[u[y,y]] returns x*y - F[u[y]], where F[u[y]] is an explicit derivative whose prettyprinted form has been shortened by declare.