[This solution is similar to MMcDara's/Sand15's "Exact solution", although I derived it independently before his was posted. I explicitly use the BetaDistribution. He uses it also, but not explicitly by name. In my (limited) experience, Statistics works better the closer one sticks to the pre-defined distributions.]

It may surprise you that your problem can be solved completely symbolically (and even using fewer lines of code than you used), such that your QStar can be expressed as a fairly simple parameterized random variable to Maple's Statistics. After that, it's trivial to compute the mean, standard deviation, etc., as symbolic expressions and to generate samples in the millions (rather than the 1000 that you asked for), all without a single loop.

Indeed, QStar can be expressed as an affine (i.e., linear plus a constant) combination of two Beta random variables. There are five parameters to the combination:

• a: the lower bound (not the sample min) of the underlying Uniform distribution
• b: the upper bound (not the sample max) of the underlying Uniform distribution
• n: the sample size (the 10 that you used)
• co: exactly as you used
• cs: exactly as you used

Here a worksheet with the analysis:
 

Download OrderStats.mw

You just need to add sin(x+t) and sin(x-t), letting x be the space parameter and t the time parameter:

plots:-animate(
   plot,
   [
      [sin(x+t), sin(x-t), sin(x+t)+sin(x-t)], x= -2*Pi..2*Pi, 
      color= [red,blue,black], thickness= [0,0,2]
   ],
   t= 0..2*Pi, size= [967,610], frames= 100
);

That's all that's needed mathematically. You can add various aesthetic doo-dads like the red dots, of course.

Here's a procedure for a linear-time^[*1], linear-memory algorithm to generate permutations:

RandPerm:= proc(n::nonnegint)
description "Fisher/Yates/Durstenfeld shuffling algorithm";
option `Reference: https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle`;
local P:= rtable([$1..n]), k, j;
   for k from n by -1 to 2 do
      j:= 1+irem(rand(), k);
      (P[k],P[j]):= (P[j],P[k])
   od;
   [seq(k, k= P)]
end proc
:

Usage:
RandPerm(9);
      [4, 6, 8, 5, 9, 3, 2, 1, 7]

A compiled-Maple version of the above procedure is possible, although I'd be skeptical about relying on a compiler-supplied version of rand. The cycle length of the generator needs to be significantly larger than n! to achieve uniformly random selection from all possible permutations. As of Maple 2018, the Maple-library version of rand has a cycle length slightly larger than 2080!, which is probably adequate for most purposes; but since 21! > 2^64, a generator based on 64-bit hardware arithmetic can't effectively shuffle a deck of playing cards, or even a half deck.

^[*1]By linear-time, I mean that the computation time is proportional to the length of the output---n in this case. Obviously, this is the best possible asymptotic time complexity.

The only practical reason I can see for doing such a computation is if it's done in modular arithmetic. Is that your overall goal? Here's an example of how it's done. The first two lines are to choose a random modulus m such as might be used in a cryptographic application.

r:= rand(2^64..2^65):
m:= (nextprime(r())-1) * (nextprime(r())-1);
          m := 868409667625403659344120587678635445760
815427339135043 &^ 5579408896058 mod m;
             28610635084817154115285869033088238569

This is returned instantaneously.

This is an example for educational purposes only. My m is still not large enough or random enough for a secure cryptographic application.

The problem with Kitonum's solution is that the root plots are discontinuous. That will almost always happen if you try to directly plot parameterized roots returned by solve (or other algebraic solution techniques). I get around that by turning each parameterized root into an IVP and solving with dsolve(..., numeric).

restart:
Digits:= trunc(evalhf(Digits)):
gm:= V -> 1/sqrt(1-V^2):
T:= w-k*V:
S:= w*V-k:
H:= unapply(
   numer(
      T*S^2*gm(V)^3*3/2+I*S^2(1+I*27/4*T*gm(V))*gm(V)^2-
      (1+I*27/4*T*gm(V))*(1+I*5*T*gm(V))*T*gm(V)
   ),
   w,V,k
):
Vlist:= [
   Record("V"= 0.8, "color"= "Black", "plot"= ()),
   Record("V"= 0.9, "color"= "Red", "plot"= ()),
   Record("V"= 0.99, "color"= "Blue", "plot"= ())
]:
R:= [solve(H(w,V,k), w, explicit)]:
d:= nops(R): #number of roots
for v in Vlist do
   S:= dsolve(
      {   #odes:
          seq(eval(diff(r[i](k),k) = diff(R[i],k), V= v:-V), i= 1..d),
          #initial conditions:
          ([seq(r[i](0), i= 1..d)] =~ evalf(eval(R, [V= v:-V, k= 0])))[]
      },
      numeric
   );
   v:-plot:= plots:-odeplot(
      S, 
      [seq([[k, Re(r[i](k))], [k, Im(r[i](k))]][], i= 1..d)], 
      k= 0..1,
      color= v:-color,
      linestyle= ['[solid,dash][]'$d]
   )
od:   
plots:-display([seq(v:-plot, v= Vlist)], size= [987,610]);

First, change all < to <=. The optimizer can't handle <. Then do:

Optimization:-NLPSolve(lhs(obj), [const[], obj], 'maximize');

Perhaps you meant sqrt(x-6) rather than sqrt(x-5)? Otherwise, there's no reason why they should "combine".

@a_simsim If you change the command copy(nodeprop) to copy(nodeprop, 'deep'), then you'll get your expected results. This is because not merely the Records themselves but also the Vectors in the Records are mutable structures, which means that a direct assignment merely creates a new pointer to the original structure rather than a new structure. The keyword deep means that copy is recursively applied to all substructures.

There is a blatant bug that can be seen in lines 5-6 of showstat(Statistics:-DataSummary): The X data themselves are being used as the weights, while the weights that you passed (as Y) are totally ignored. It's exactly as if it were behaving correctly and you had specified DataSummary(X, weights= X).

The command is GraphTheory:-SetVertexPositions.

X:= Matrix(9, 9, [[0, 0, 1, 0, 1, 1, 1, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0]]):

#Coordinates
P:= [[40, 40], [22, 22], [36, 26], [21, 45], [45, 35], [55, 20], [55, 45], [26, 59], [55, 65]]:

GT := GraphTheory;
G := GT:-Graph(X);
GT:-SetVertexPositions(G, P);
GT:-DrawGraph(G);

It is interesting that the position information is stored in the graph itself rather than simply in the plot. Also, note that SetVertexPositions modifies its first argument, the graph, rather than giving a return value.

You need to use all vertices, in order, in a list; so, tweaking is best done by alternating calls to GetVertexPositions and SetVertexPositions. Several examples are shown at help page ?GraphTheory,GetVertexPositions.

Is this what you're looking for?

restart:
A:= <
   1.364, 0.64765768e-1  ; 
   2.05,  -.182176113    ; 
   2.664, -0.13914542e-1 ; 
   2.728, 0.2193938e-1   ; 
   4.092, -0.18349139e-1 ; 
   4.1,   -.312968801    ; 
   5.328, -0.1819837e-2  ; 
   5.456, -.28840961     ; 
   6.15,  -.57076866     ; 
   7.992, .175022254
>: 
F:= unapply(CurveFitting:-LeastSquares(A, x), x):
P:= plot(
   [F, A], 0..8, color = [red, blue], style = [line, point],
   legend= ["Метод наименьших квадратов", "Экспериментальные данные"], 
   legendstyle= [font = ["Roman", 15]],
   labels= [typeset(d, ", ", `мм`), typeset(ln(`I__ `/I__0))], labelfont = ["Roman", 15],
   axesfont= ["ROMAN", "ROMAN", 15], linestyle = [solid], symbolsize = 20, 
   title= "Определение линейного коэффициента поглощения", titlefont = [Roman, bold, 20]
);
X:= A[..,1]:  Y:= A[..,2]:  `Y^`:= F~(X):
plots:-display(
   [
      P,
      Statistics:-ErrorPlot(
         `Y^`, xcoords= X, yerrors= abs~(Y-`Y^`), markers= false, thickness= 0
      )
   ],
   axes= normal, gridlines= false
);

Rather than putting an upper bound on k in the for loop, you could put an upper limit on the time your willing to spend on the computation. Here's my version:

Wrapped_prime:= proc(p::prime, {base::posint:= 10, filter::anything:= (()-> true)})
global lastK;
local pow:= base^trunc(ln(p)/ln(base)), b2:= base^2, k, m:= p;
   lastK:= 0;
   for k do
      #Recursive computation of repunit-wrapped number:
      pow:= b2*pow;  m:= pow + base*m + 1;
      if filter(k) then 
         if isprime(m) then return k fi;
         lastK:= k
      fi
   od
end proc
:     

Since repunits exist in any base, not just base 10 (Mersenne numbers being the most-famous example), the procedure will work in any base, defaulting to 10. The keyword parameter filter lets you specify a filtering of the k values, for your example, it's filter= (k-> irem(k,6) = 0). The  global lastK lets you check what was last value of k for which the wrapped number was verified composite by isprime.

The first argument to timelimit is the number of seconds of CPU time to expend before giving up.

timelimit(5, Wrapped_prime(397, filter= (k-> irem(k,6)=0)));
Error, (in isprime) time expired
lastK;
                              1656
timelimit(5, Wrapped_prime(397, base= 5));
                               2
timelimit(5, Wrapped_prime(397, base= 2));
                               24

Your system of equations contains numerous terms of the form

sum(FyjR, j = 1 .. m),

which makes it seem as if you believe that FyjR can be indexed by j; it can't be. That sum evaluates to simply m*FyjR.

Maple's mod is built to do computations in higher algebraic structures (such as rings of polynomials and matrices) built on top of the ring of integers modulo m. Thus, unevaluated variables such as your n are treated as transcendentals (or "dummy variables") in a polynomial, so n mod 2 is truly equal to n. In your case, you want n to represent an integer (which presumably will be supplied later), not a transcendental. In that case you should use irem(n,2) instead of mod(n,2).

While I strongly recommend that you use irem instead mod for this situation, Tom Leslie shows that you can use mod(n,x) if both n and x are unevaluated and integer values will be supplied for both simultaneously, or at least n is supplied before x.

Okay, you're running the command-line version of Maple interactively. Let's suppose that the file that you created with Vim is named "C:/dir1/dir2/MyCode.mpl". Then, at the ">" cursor/prompt, enter

read "C:/dir1/dir2/MyCode.mpl";