@testht06 

Firstly, I've streamlined the procedures Lehmer and LehmerInv by using a remember table for factorials and using 1-relative indexing for the coefficients c[i]. Here are the new procedures:

#This is simply a remember table for factorials:
Fact:= proc(n::nonnegint) option remember; n*thisproc(n-1) end proc:  Fact(0):= 1
:
Lehmer:= proc(S::And(list(nonnegint), satisfies(S-> max(S)=nops({S[]})-1)))
description 
   "Rank (0..nops(S)!-1) of a permutation of [$0..nops(S)-1] using the"
   " canonical Lehmer order."
;
local i, s, n:= nops(S);
   add(Fact(n-i)*add(`if`(s < S[i], 1, 0), s= S[i+1..]), i= 1..n-1)
end proc
:
LehmerInv:= proc(s::nonnegint, n::And(posint, satisfies(n-> s < Fact(n))))
description "Returns a permutation of [$0..n-1] from its Lehmer code.";
option `Reference: https://en.wikipedia.org/wiki/Lehmer_code`;
local i, c:= s, C:= rtable(1..n, i-> iquo(c, Fact(n-i), 'c'));         
   for i from n by -1 to 2 do
      c:= C[i-1]; 
      C[i..]:= C[i..] + rtable(i..n, j-> `if`(C[j] < c, 0, 1))
   od;
   [seq(c, c= C)]
end proc
:

Now, onto your sequence computation. Note that the composition of permutations is not commutative, so your "x" operator is ambiguous. Below, I use both possible interpretations. Also, n needs to be a parameter in addition to c and x0.

MySeq:= proc(c::realcons, x0::nonnegint, n::And(posint, satisfies(n-> x0 < Fact(n))))
local
   N:= Fact(n)-1, P:= Pi/2, i, C:= (c+1)/(N+2), l:= Lehmer, L:= LehmerInv,
   MySeq:= proc(i) 
   option remember; 
      (x-> l(L(x,n) &x L(floor(N*sin(P*(x/N+C))), n)))(thisproc(i-1))
   end proc
;
   MySeq(0):= x0;
   seq(MySeq(i), i= 0..N)
end proc
:
`&x`:= (P,Q)-> P[Q+~1]: #Compose permutations left first.
MySeq(4, 5, 4);
5, 16, 6, 13, 7, 23, 2, 16, 6, 13, 7, 23, 2, 16, 6, 13, 7, 23, 2, 16, 6, 13, 7, 23

`&x`:= (P,Q)-> Q[P+~1]: #Compose right first.
MySeq(4, 5, 4);
5, 16, 1, 9, 0, 7, 23, 2, 7, 23, 2, 7, 23, 2, 7, 23, 2, 7, 23, 2, 7, 23, 2, 7

So, the first sequence is attracted to a 6-orbit, and the second to a 3-orbit.

@testht06 Change the second-to-last line from [seq(C)] to [seq(c, c= C)].

@testht06 

For example,

LehmerInv(13, 4);
      [2, 0, 3, 1]

Note that the length of the permutation, 4 in this case, can't be determined solely from its Lehmer code, 13. Compare with

LehmerInv(13, 5);
      [0, 3, 1, 4, 2]

Let me know if it's working for you. There could be some version issue that I missed and that I have no way to test for.

Will the matrix dimension always be 7x7?

Here they are as straightforward procedures direct from the algorithms, no packages needed. The primary procedure, Lehmer, I coded directly from your Question. For the inverse, I used an algorithm given in the Wikipedia article "Lehmer code". (I will look up Knuth's algorithm when I get out of bed :-); the double loop is a little scary.)

Lehmer:= proc(S::And(list(nonnegint), satisfies(S-> max(S)=nops({S[]})-1)))
description 
   "Rank (0..nops(S)!-1) of a permutation of [$0..nops(S)-1] using the"
   " canonical Lehmer order."
;
local i, n:= nops(S);
   add((n-i)!*nops(select(s-> s < S[i], S[i+1..])), i= 1..n-1)
end proc:

LehmerInv:= proc(s::nonnegint, n::And(posint, satisfies(n-> s < n!)))
description "Returns a permutation of [$0..n-1] from its Lehmer code.";
option `Reference: https://en.wikipedia.org/wiki/Lehmer_code`;
local i, q:= s, C:= rtable(0..n-1, i-> iquo(q, (n-i-1)!, 'q'));
   for i from n-2 by -1 to 0 do
      q:= C[i]; 
      C[i+1..]:= C[i+1..] + rtable(i+1..n-1, j-> `if`(C[j] >= q, 1, 0))
   od;
   [seq(C)]
end proc:

Note that the inverse procedure requires a second argument, n, the permutation's length.

@testht06 The Iterator combinatorics package was introduced in Maple 2016. If you're stuck with Maple 15, I'll need to write something compatible with that.

@David Sycamore My comment was about the commands sqrfree and isqrfree. Yes, capitalization and spelling matter!

@testht06 Lehmer([0,2,1,3]) should work and return 2. So, you get nothing, not even an error message? That's weird. What Maple version are you using?

What do you want it to return for x*y?

What does "2" mean in your (dx/dt)2?

@Carl Love I also attempted a quadratic residue analysis to prove that there were no other solutions. This, however, failed: I was easily able to find at least one quadratic residue (other than f(1), of course) for every modulus up to 2^17 = 131072. Here's my code:

F:= rsolve({f(n) = 4*f(n-1) - 2*f(n-2), f(1)=1, f(2)=28}, f(n), makeproc):
upperN:= 2^13: upperM:= 2^17:
FL:= F~([$2..upperN]):
for m from 2 to upperM do
   if andmap(x-> numtheory:-quadres(x,m) = -1, FL) then print(m); break fi
od:

@acer This is a great Answer, spectacular in several senses of that word!

A Post that I wrote (which you made significant contributions to) within the past year delves significantly into dropped-shadow contour plots (such as yours above) generated from parameterized  dsolve(..., numeric, method= bvp) output: "Numerically solving BVPs that have many parameters". (I know that you're already aware of this; I mention it for the benefit of other readers.)

@Jjjones98 The code above is a procedure definition; it is not intended to have any output. To get output, you need input, specifically the matrix size; you need to enter something like UMG(5).

@mmcdara 

In my opinion, those extra characters that may appear at the beginnings of continuation lines within an execution group that look the same as the prompt characters are just garbage, bugs, that only occur in some versions of Maple. They serve no useful purpose, are a total nuisance to remove, and I don't consider them to be prompts because a prompt is what your cursor sits to the right of while awaiting input. In other words, I'm saying that there's no such thing as a "passive" prompt; there are only active prompts and garbage characters.

However, I understand your point that my parenthetical comment could be misunderstood because those garbage characters do look similar to prompts; however, they can be distinguished from true prompts by paying attention to the fine, black bracket lines in the left margin.

@David Sycamore Kitonum's code checks the square-free property by checking that all exponents in the prime factorization are 1. The commands sqrfree and isqrfree are not for checking the condition. Rather, they are for producing a square-free factorization. This is a bit more elaborate, and it'd be wasteful to use these for the purpose of checking the condition.