Please upload (at the very least) the code, the call that produces the error, and the exact error message. You can either upload a worksheet as an attached file, or just upload plaintext. It doesn't matter if it's in some other language; we can still figure it out.

@tarik_mohamadi 

Is the case R=0 meaningful to you? It certainly wouldn't be meaningful to me. This case is trivial and degenerate, and it needs special handling if you want to include it. Otherwise, just remove this case from the loop.

@tarik_mohamadi I need to see your complete code to provide further help. At the least, I need to see the value of eqns.

Generally, if it's possible to do something with Threads, then doing so is significantly faster and uses significantly less memory than using Grid. So, if some easy modification can make the code threadsafe, it may be worth pursuing that.

It's possible for code to be threadsafe even though it uses procedures not listed at ?index,threadsafe. That page lists procedures that are guaranteed to be threadsafe. The issue is whether the parallelized part of the code can make contradictory assignments to variables at a higher lexical level (including enviroment variables such as Digits). A common case is code running under evalhf. I believe that such is always threadsafe (and, if not, it certainly is threadsafe in the vast majority of cases). So, if that evalhf code uses sin, for example, it's still threadsafe even though sin, in general, is not.

@Joe Riel According to the relevant help page, combinat:-rankperm and combinat:-unrankperm were introduced in Maple 16; so, too late for this OP. And, I'm content with the code that I've written, but if you see the possibility of improvement, please let me know. I'm aware that with Maple 2019 syntax, the line

C[i..]:= C[i..] + rtable(i..n, j-> `if`(C[j] < c, 0, 1))

can be improved to 

C[i..]+= rtable(i..n, j-> `if`(C[j] < c, 0, 1))

@testht06 

My composition operator is the procedure &x that I already gave above:

[3,5,7,0,2,1,4,6] &x [2,4,0,5,3,1,6,7];
      [5, 1, 7, 2, 0, 4, 3, 6]

Maple's sophisticated list indexing makes the code of this procedure so short that you may have missed it above:

`&x`:= (P,Q)-> Q[P+~1]:

@testht06 

Okay, your composition of permutations is what I called "compose right first", and that's how I left it in the following.

Once the sequence's orbit has been found, it's redundant to continue generating it. So, I've replaced MySeq (procedure above) with FindOrbit, which stops as soon as the orbit is found. Using c=4, x0=5, n=8 the orbit length is 177 (and is returned in about 1/20 of a second on my computer).

FindOrbit:= proc(c::realcons, x0::nonnegint, n::And(posint, satisfies(n-> x0 < Fact(n))))
local
   N:= Fact(n)-1, C:= (c+1)/(N+2), P:= Pi/2, l:= Lehmer, L:= LehmerInv, S, O, i, j, x,
   It:= x-> l(L(x,n) &x L(floor(N*sin(P*(x/N+C))), n))
;
   (x, S[0], O[x0]):= (x0, x0, 0);
   for i do
      x:= It(x);
      if O[x]::nonnegint then return i-O[x], [seq(S[j], j=  0..i-1)] fi;
      (S[i], O[x]):= (x, i)
   od         
end proc
:
`&x`:= (P,Q)-> Q[P+~1]: #Compose right first.
:
(Len, Seq):= CodeTools:-Usage(FindOrbit(4, 5, 8));

Len, Seq := 177, [5, 13, 30, 24, 105, 196, 590, 766, 3782, 8279, 20411, 13524, 
13496, 9576, 39969, 6110, 37175, 1470, 1749, 4636, 9958, 3331, 10039, 1167, 3499,
7082, 7923, 24720, 7913, 22906, 3029, 1751, 4831, 5059, 21591, 9385, 38428, 17012, 
26884, 27454, 34552, 20239, 39443, 996, 1093, 2568, 182, 475, 1243, 4224, 8551, 
29204, 17289, 15468, 12379, 36538, 18721, 37976, 879, 4779, 9236, 30551, 19993,
519, 1305, 4491, 8508, 27959, 14216, 15316, 38489, 6888, 557, 1083, 2817, 1752, 
4865, 6503, 2285, 2132, 4539, 8555, 28483, 17795, 24620, 23479, 23266, 37753, 
26326, 13309, 26102, 13299, 25398, 19246, 4983, 5500, 25942, 39541, 11698, 30111,
15606, 33425, 12072, 10824, 36352, 4119, 6599, 3017, 358, 655, 1094, 2585, 462, 
1014, 1310, 4896, 5760, 32619, 25650, 39754, 7045, 5731, 31460, 11862, 39758, 7041,
5747, 31434, 6641, 285, 698, 1434, 207, 562, 1049, 1381, 455, 1295, 5012, 9173, 
31495, 11530, 9649, 39249, 20510, 8557, 28641, 25632, 9115, 30301, 9526, 37674, 
16403, 19310, 32361, 10856, 36255, 4073, 7755, 16182, 9221, 31378, 6559, 4482,
8480, 27918, 6917, 2232, 427, 16, 28, 6, 23, 32, 42, 8, 7]

The returned values are the orbit length followed by the entire sequence upto (but not including) the first repeated element. It is just a coincidence in this case that the entire sequence is the orbit.

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