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

@acer So, multiplying both sides by the denominator is all that's needed to get rsolve to explicitly solve it!

@acer Thank you. So, it seems that the inclusion of rtable_eval incurs a slight performance penalty (see test code below), while not doing anything useful. 

Here's my motivation for these Questions. Consider the following simple and pedagogically super-sweet little procedure for computing the GCD and Bezout coefficients of any two integers using a matrix-based method:

GCDwithBezout:= proc(a::integer, b::integer)
local q, r0:= a, r1:= b, S:= <0,1; 1,0>, Q:= <0,1; 1,-q>;
     while r1 <> 0 do
         (r0, r1):= (r1, irem(r0, r1, 'q'));
         S.= eval(Q)
     od;
     sign(r0)*[r0, seq(S[..,1])]
end proc:

I realize that this can be done much faster with igcdex, but I am considering generalizations which can't be.

I've noticed that I can get a slight performance boost by declaring S with datatype= integer. This makes sense, and it's done for all the examples below.

I'm interested in finding the fastest way to do the product update S.= eval(Q). So far, that itself is the fastest way that I've found. Here are 5 alternatives, with the very suprising result that inplace operation is by far the slowest method:

GCDwithBezout:= proc(a::integer, b::integer)
local 
     q, r0:= a, r1:= b, Q:= <0,1; 1,-q>,
     S:= rtable((1..2)$2, [[1,0],[0,1]], 'datatype'= integer, 'subtype'= Matrix)
;
     while r1 <> 0 do
         (r0, r1):= (r1, irem(r0, r1, 'q'));
         #(* Alternative 1:
              S.= eval(Q)
         #*) 
         (* Alternative 2: Same time as #1:
              S:= S.eval(Q)
         *)
         (* Alternative 3: Takes a bit longer than #1 or #2:
              S.= rtable_eval(eval(Q))
         *)
         (* Alternative 4: Takes almost twice as long as #1 or #2:
              S:= LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(S, eval(Q))
         *)
         (* Alternative 5: Takes 7-8 times (!!!) longer than #1 or #2:  
              LinearAlgebra:-LA_Main:-MatrixMatrixMultiply(S, eval(Q), 'inplace')
         *)
     od;
     sign(r0)*[r0, seq(S[..,1])]
end proc:

And here is a test suite for it, or feel free to construct your own:

#Code to test it:
#Fast, semi-builtin code, for comparison:
Via_igcdex:= proc(a,b) local s, t, g:= igcdex(a,b,s,t); [g, s, t] end proc:

#Procs being tested:
Procs:= [GCDwithBezout, Via_igcdex]:

#Generate random cases to test:
#(words per integer; number of test cases; random seed):
(W, N, s):= (2, 2^12, 1):
randomize(s):
R:= rand(map(`*`, -1..1, 2^(W*kernelopts(wordsize)-1)-1)):
(A, B):= ('['R()' $ N]' $ 2):

#efficiency test:
R:= CodeTools:-Usage~(map(''curry(printf, "\n%a: "), zip'', Procs, A, B)):

GCDwithBezout: memory used=78.44MiB, alloc change=0 bytes, 
cpu time=671.00ms, real time=626.00ms, gc time=140.62ms

Via_igcdex: memory used=7.03MiB, alloc change=0 bytes,
 cpu time=47.00ms, real time=59.00ms, gc time=0ns

nops({R[]}); #accuracy test (1 means all Procs agree).
                               1

All of the methods give accurate results; efficiency is the only issue here.