If you're using Maple 2019 or later, you can do this:

A:= Array([4]);
a:= A[-1]: to 3 do A,= ++a od;

The ,= is the new "append" operator, which is especially useful for efficiently appending to Arrays. 

The ++ is a new operator that does exactly what it does in C++: ++a adds 1 to a, stores the result in a, and returns the updated value.

The A[-1] gets the last element of A (which is also the first element in this case). This is not new syntax.

You can use seq like this:

A:= Array([4]):
n:= 3:
u:= upperbound(A): 
a:= A[-1]: A(u+1..u+n):= <seq(a+i, i= 1..n)>;

To me, that's nowhere near as clear as the new syntax.

All of the above assumes that you want to append to the end of the Array. If instead you want to insert new elements somewhere other than the end, that's still possible to do, but the commands required are quite different. Let me know if that's what you want. 
 

It seems that the piecewise expressions returned by PDF confuse the numeric integrator's ability to make an accuracy estimate. The simplify command applied to the integral (in Int or int ​​​​​​form) will select the relevant branch of the piecewise, and the numeric integration will finish.

Your equation has its variable, q, under the square root. As you may recall from high-school algebra, solving such an equation involves first converting it to an equivalent polynomial equation. Some of the polynomial's solutions may not work in the original equation because the square root indicates only the positive branch. Because your equation also has the parameters c, m, and s, there is no single expression for q that will work for all values of the parameters. That's why back-substituting and simplifying doesn't give you 0.

There is, however, an elaborate process by which you can classify these conditional solutions and know what are the conditions that the parameters must satisfy for each to be valid. It uses the parametric option to the solve command. I show it in the attached worksheet.

Download SolveParametric.mw

Here are the highlights of the worksheet:

restart:
p:= (1-u)*(u-m)-q*(u+c)+u*(1+m-2*u)-s:
u:= (1+m-q-sqrt((1+m-q)^2-4*(m+c*q)))/2:
p:= simplify(p);

The command evala(Norm(p)) finds the minimal polynomial equivalent to p. In other words, it removes the radical. 

pp:= evala(Norm(p));
         2  2      2            2                  2          
pp := 4 c  q  - c m  q + 2 c m q  + 6 c m q + 2 c q  - 4 c q s

      3      2      2                  2              2        
   - m  + 2 m  q + m  s - 3 m q s + 2 q  s - c q + 2 m  + 2 m q

                      2        
   - 2 m s - 3 q s + s  - m + s

I will assume that you're interested only in real solutions. Thus, the expression under the radical must be nonnegative. The parametric option to solve finds conditions on the parameters (c, m, and s in this case) that make certain solution branches correct. It only works on polynomials.

The following command takes about 30 seconds.

pwsol:= solve({pp, q^2 + (-4*c-2*m-2)*q + (m-1)^2 >= 0}, q, parametric, real);
nops(pwsol);
                              1731

That means that there are 866 cases! I list them in the worksheet.

My preferred alternative to Acer's 'SF' is :-SF. The :- prefix unmistakably marks it as global (to a Maple-savvy reader). Preben's global SF also unmistakably marks it as global, of course, but I prefer the :- because it's targetted to the specific instance of the variable and not other possible occurences in the same procedure.

':-SF' also works.

First, data and d need to be Vectors, not Matrixs!

(`+`@seq/numelems)~(
    map2(
        `?[]`, 
           d, 
        `[]`~(remove(`=`, [ListTools:-Split(j-> data[j]>0.01, [$1..numelems(data)])], []))
    )
); 
             [22/3, 51/2]

Here's one way:

{seq(seq(seq(seq((a*x0+b*x1)*(c*y0+d*y1), a= {0,1}), b= {0,1}), c= {0,1}), d= {0,1})};

Of course the result is 0 for the reasons given by the others. Here's a nice way to get Maple to say that directly:

f:= z-> z/abs(z)^2:  r:= t-> exp(I*t):
IntegrationTools:-Change(Int(f(z), z= r(-%Pi)..r(%Pi)), z= r(t));
value(%);

Notice how % is used to make Pi inert, the inertness being removed by value.
 

Many algorithms use some sort of randomization. For example, the initial values of the decision variables could be chosen at random. For several purposes---for example, debugging or documentation---it's useful to be able to rerun code using the same random values. To do this, use the randomize(...) command. For example:

randomize(42); 
DirectSearch:-Search(...);

The 42 could be any positive integer. As long as it's the same number on different runs, the random results should be the same, even if it's run on a different computer. If some program doesn't respect this (for example, by using an ad hoc random generator), it should be considered a bug. I haven't checked whether DirectSearch respects this.

If you want to vary the results but still have repeatability, do 

key:= randomize(); #no arguments
             key := 1436586243734
DirectSearch:-Search(...);

Now to repeat using the same random numbers, do

randomize(1436586243734): 
DirectSearch:-Search(...);

If you've not done a restart, the long number can be replaced by key.

Any such repetitive code can be and should be simplified. Your whole code above can be replaced by

V:= [E,L,P,M]:
Vf:= (op@index~)~([V,V], [1,2]);
dV:= [d||~((V||~1)[], (V||~2)[])]; 
Data:= {a= 1, b= 2, c= 3}: #Put your actual data here.
Sistem:= diff~(Vf(t),t) =~ subs(Vf=~ Vf(t), Data, dV);

The direct cause of your error is that Data must be an equation or a list or set of equations (as shown above).
 

Like this:

eq:= diff(z(x,t),t) + diff(R(x,t),x)*L(x,t) + 
    A*diff(L(x,t),x) + diff(k(x,t),x,x) - diff(k(x,t),t,t) = 
    0;
PDEtools:-dchange(
    eliminate({X= x-alpha*t, T=beta*t}, {x,t})[1],
    eq, [X,T], params= {alpha, beta}
);
 

Your procedure only makes sense for integer n, but plot, by default, passes it real values in the interval 0..10.

One possible correction is 

plot(b@trunc, 0..10)

The trunc will convert the numbers passed to b by plot into integers. The @ is the function composition operator.

There are numerous other ways to correct this.

You could use FileTools:-ModificationTime.

The one-line procedure AssignToClusters does both jobs.

To create a discrete distribution with support other than positive integers, use Empirical rather than ProbabilityTable. For example,

pr:= [1/2, 1/4, 1/8, 1/16, 1/32, 1/32]:
X:= Statistics:-RandomVariable(
    Empirical([0,1,2,3,4,5],
    probabilities= pr)
):
S:= Statistics:-Sample(X, 1000):
plots:-display(
    Statistics:-Histogram(S, discrete, thickness= 5),
    plot([seq([h,pr[h+1]], h= 0..5)]),
    axes= normal
);

 

Do

solve(ineq1, chi, parametric);