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

Your definition of ff in the worksheet refers to a Matrix ffd, which doesn't seem to exist in the worksheet. Thus, the worksheet can't be properly executed.

@tomleslie Thanks. Perhaps I was too quick in saying "ignore the final paragraph": It's essential to my purpose that the rtable first have a symbolic value (the x in this case), and only then that that variable be assigned a value. My other Question posted yesterday, "Unexpected difference...", shows the reason, although it's still an over-simplification.

I want to compute the product of a large number of small matrices of arbitrarily large integers. The matrices differ only in one element. Of course, I can easily figure out how to do this; what I'm interested in is doing it in the most efficient way possible. My idea was to put x in the matrix's variable position, and in a loop assign the values to x and accumulate the product. There should be a way to do this and only create a total of three matrices: one for the matrix with the symbolic entry, one for its instantiations--which'll get overwritten rather than created anew, and one for the accumulated product.

@acer Thanks. What do you think of the sentence immediately after what you quoted, which says "Applying eval() to an rtable returns immediately without changing the given rtable."? But your example shows that it must do something.

@Kitonum Thanks. Only the first of your three evals is needed; I think that the others do nothing at all.

Of course, my Question still stands: Why the difference, etc.?

@Kitonum Thanks. I did come up with that as a solution to my actual problem.

I found something else, almost certainly a bug, that is the root cause of my problem. So, please ignore the final paragraph of the Question.

@TeyhaNeedHelp I'm not sure what you mean by "display the value in numeric". Do you mean that you want to see 6 on the axis rather than the infinity symbol? If so, then simply remove the xtickmarks option (the final option) from my command.

I'm also not sure what you mean by "just use dsolve". The dsolve command does not directly produce plots; its output needs to be sent to a plotting command.