There is a simple bug in `evalf/doublefactorial`(z) only for z=1. The bug is in line 5 of showstat(`evalf/doublefactorial`): The 1.+z should be 1.. (And I wonder why whoever wrote this thought that z=1 needed to be treated as a special case. It seems to me that only z=0 needs special treatment.)

You can work around the bug by setting 

`evalf/doublefactorial`(1.):= 1.:

at the start. This is called a remember-table assignment.

This bug only affects numerical evaluation (via using evalf or floats). Symbolics aren't affected.

These three things that you believe are either not true or not valuable, and if you did accomplish them (which I could easily tell you how to do), you'd find out after that they aren't as good as the better ways that I show below:

1. "But since the list contains the values as an expression, for example DQV[0, 0, 0, 0] = 0. In such cases, I cannot perform mathematical operations for the solution points results."
2. "if I can delete the DQV[0,0,0,0] from the expression mentioned above, I'll end with the numeric values list, which can then be converted into an array and make things easier for me."
3. "I'm struggling to delete the DQV part of the expression in one go (manually deletion for each expression, it's very time consuming to delete for 1500 variables)."

The first isn't true because you can easily use eval(expr, S) (where expr is any expression containing any subset of the variables and S is defined below) to numerically evaluate any mathematical expression of the variables. An example of using eval is given in example 3 below.

The second isn't valuable because it can easily be done without deleting anything, as shown in examples 2 and 4 below.

The third isn't valuable. In other words, there's no need to manually edit any output. I can't imagine anything useful that could be done with an LPSolve solution that couldn't be done with some simple commands (usually just 1 line) and no manual editing.

#Extract the list of variable = value equations:
S:= A:-Results("solutionpoint"):

#Example 1: Find all variables that are equal to 0:
ZeroVs:= lhs~(select(e-> rhs(e)=0, S));

#Example 2: Create a table (essentially an array) of all DQV variables
#(You shouldn't name it DQV, but Dqv is fine):
Dqv:= table(((op@lhs)=rhs)~(select(e-> type(lhs(e), specindex(DQV)), S))):

Dqv[1,1,1,0]; #etc.

#Example 3: Let C be your set of constraints. The following finds the subset that are
#inequality constraints that are at their limit at the current solution:
ActiveCs:= select(c-> c::`<=` and eval((lhs-rhs)(c), S) = 0, C);

#Example 4: Like example 2, but create a table of tables of all variables:
All:= table():
for e in S do
    if lhs(e)::indexed then All[op(0,lhs(e)][op(lhs(e))]:= rhs(e)
    else All[lhs(e)]:= rhs(e)
    fi
od:

All[DQV][1,1,1,0]; #etc.

I'd be happy to give more such examples if you say more precisely (mathematically) what you want. And you should ask here if you ever feel a need to manually edit any Maple output to produce input for other Maple commands.

Change the for line to

for ind, val in eval(F_vRk) do

Tables, procedures, and some modules have a property called last name evaluation. See ?last_name_eval.

Download RSA.mw

There are several ways to construct the polynomials. I think the easiest is to use the command quo, which divides two polynomials, returning a polynomial quotient (and, optionally, a polynomial remainder, but it's always 0 in this case).

restart:

unimodal:= proc(L::list(numeric))
local n:= nops(L), k;
    for k to n-1 while L[k] <= L[k+1] do od;
    for k from k to n-1 while L[k] >= L[k+1] do od;
    evalb(k=n) #i.e., Did we get to the end?
end proc
:
qbinomial:= (n::nonnegint, k::nonnegint, q::name)-> 
    local i; 
    quo(mul~(q^~[n-k+i, i] -~ 1, i= 1..k)[], q)
:

qbinomial(6,3,q);
             

A good way to count the things that satisfy a condition is to add 1 everytime you find one, like this:

PT:= PolynomialTools:
add(
    add(
        `if`(unimodal(PT:-CoefficientList(qbinomial(n,k,q),q)), 0, 1),
        k= 0..n
    ), 
    n= 0..10
);
                               0

Given acer's Comment, I hope that I don't get into trouble for Answering this. Nonetheless, I agree with what he said. Anyway, the Answer is trivial, so here goes:

Your line 

_Envdiffopdomain:= [Dt,t]:

specifies that your special^[*1] differential operator is Dt. But your expression tries to use D, not Dt, as the differential operator. That's all there is to it.

[*1] By special, I mean that Dt is the differential operator that diffop2de recognizes as the one that it's supposed to treat specially in the context of diffop2de. 

It's shocking to me that your instructor proposes an extremely inefficient O(n^2) solution when it's obvious even to a child that there's an easy-to-understand O(n) solution. (Of course, the child wouldn't think of it in those words; they'd just know it intuitively.) In other words, it's obvious that unimodality can be determined by making a single pass through the list, and when the result is false, that can usually be determined before reaching the end. Anyway, here's my O(n) procedure:

unimodal:= proc(L::list(realcons))
local n:= nops(L), k;
    for k to n-1 while is(L[k] <= L[k+1]) do od;
    for k from k to n-1 while is(L[k] >= L[k+1]) do od;
    evalb(k=n) #i.e., Did we get to the end?
end proc
:

Hint 2: None of that given information about phi (also called the totient) is needed. There's a very easy way to solve this just knowing that p and q are integers greater than 1 and that
             p*q = 10^100 + 598*10^50 + 67497

I have more hints after you think about that one for a while. 

Matrices can be entered with the rows separated by semicolons:

M:= <
    3626,  -189,   -513,   -630,   -945,   -1323  ;
    14544, -1392,  -2784,  -4415,  -6960,  -9744  ;
    12992, -2016,  -4032,  -6720,  -10215, -14112 ;
    81408, -19200, -38400, -64000, -96000, -133455;
    26,    -9,     -18,    -30,    -45,    -63    ;
    6,     -3,     -6,     -10,    -15,    -21    #
>;

This works in both 1D and 2D Input. @Scot Gould showed a similar method using the ASCII vertical bar as the separator; however, with that method, one is essentially entering the transpose of the desired matrix. By using the semicolon separator, you get a matrix with exact visual correspondence to your input.

Perhaps this will be sufficient for you. I changed the comma to the American decimal point. Just type in these commands, the lines in red. The black part is the output, but will be blue and gray in Maple.

S:= Sum(S__k^2*v__k, k= 1..n) = 71.70825:
print(S);
   

The styling isn't exactly the same as what you posted. I'm just starting with the easiest option in case that's good enough for you.

For a linear, homogenous system of ODEs such as yours, any Runge-Kutta method can be done as a single matrix-vector multiplication for each solution point: the previous solution point times the matrix. If the system is also autonomous (such as yours), the matrix doesn't change; if it's not autonomous, the matrix can be expressed as a matrix function of the independent variable. The amount of computational effort needed to construct these matrices is trivial. So, if the number of solution points is large, all R-K methods applied to systems such as yours require essentially the same amount of effort regardless of the order of the method.

You said that you didn't see a reduction in error when going from RK2 to RK3 to RK4. My results below contradict that; I get a substantial error reduction at each higher order.

restart:
Digits:= 15:
LA:= LinearAlgebra:

#Setup ODEs. This is your system entered in a neater form.
As:= [seq](A[k], k= 2..7)(t): #dependent functions
nsys:= nops(As):
M:= <
    3626,  -189,   -513,   -630,   -945,   -1323  ;
    14544, -1392,  -2784,  -4415,  -6960,  -9744  ;
    12992, -2016,  -4032,  -6720,  -10215, -14112 ;
    81408, -19200, -38400, -64000, -96000, -133455;
    26,    -9,     -18,    -30,    -45,    -63    ;
    6,     -3,     -6,     -10,    -15,    -21    #
>:
pi:= [seq](Pi^k, k= 0..nsys-1):
sys:= diff~(As,t) =~ 
    [seq](x, x= M.<As*~pi>)*~[-28, 28, -70, 14, -28672, 32768]/~pi/~(315*Pi^2):

#initial conditions:
ICs:= <
    -0.001210685373, -0.1636380032,   -0.003838287851,
    0.01100619795,   -0.001005637627, -0.00002775982849
>:
(t0,tf):= (0,1): #solution interval
npts:= 21: #evaluation points (including ICs at t0) 
for k to 4 do #create matrices to hold solutions for RK1, ..., RK4
    resMat[k]:= Matrix((npts,nsys), datatype= hfloat):
    resMat[k][1]:= ICs:
od:
h:= evalf((tf-t0)/(npts-1)): #stepsize

#Extract coefficient matrix of ODE system and transpose it: 
f:= Matrix(evalf[18](LA:-GenerateMatrix(rhs~(sys), As)[1])^%T, datatype= hfloat):

#Construct matrices for RK1, ..., RK4:
Id:= Matrix(nsys, shape= identity):
RK[1]:= Id+h*f:          #aka forward Euler's method 
RK[2]:= (Id+RK[1]^2)/2:  #aka Heun's method          
RK[3]:= (Id+(Id+RK[2]).RK[1])/3:                   
RK[4]:= (3*Id + (2*(Id+RK[3])+RK[1]).RK[1])/8:      

#Solve by each method:
for k to 4 do
    for j from 2 to npts do resMat[k][j]:= resMat[k][j-1].RK[k] od 
od:
     
#Get matrix-form solution from high-accuracy standard Maple solver:
sol:= dsolve(
    {sys[], seq(x, x= (eval(As, t= t0)=~ICs))}, numeric, abserr= 1e-13,
    output= Array(1..npts, k-> t0+(k-1)*h)
)[2,1][.., 2..]
:
#Compute sum-squared error (all 6 A(t) functions) for each method:
for k to 4 do `Sum-squared error RK`[k] = add(x^2, x= sol-resMat[k]) od;
        Sum-squared error RK[1] = 0.0000302404426180997
                                                     -8
        Sum-squared error RK[2] = 9.36429863900565 10  
                                                     -9
        Sum-squared error RK[3] = 6.24726072793177 10  
                                                     -10
        Sum-squared error RK[4] = 2.24644854174923 10   

Maple has at least three packages for constructing combinatorial objects: combinat, Iterator, and combstruct. Additonally, many combinatorial objects can be created by nesting (or folding) seq commands to an arbitrary depth (by using foldl). This method tends to be extremely fast. The code below is 4-5 times faster than Tom's combinat method and my polynomial method (both of which take roughly the same time) when there are thousands of sublists or more.

AllS:= proc(L::list, n::posint)
local j, k, v:= nops(L), K:= [seq](k[j], j= 1..v);
    [value(
        {foldl}(%seq, `$`~(L,K), seq(k[j]= 0..n-`+`(K[j+1..][]), j= 1..v))
    )[2..][]]
end proc
:
L__s:= CodeTools:-Usage(AllS([x,y,z], 4));
memory used=17.08KiB, alloc change=0 bytes, cpu time=0ns, real time=0ns, gc time=0ns
L__s := [[x], [y], [z], [x, x], [x, y], [x, z], [y, y], [y, z], 
  [z, z], [x, x, x], [x, x, y], [x, x, z], [x, y, y], [x, y, z], 
  [x, z, z], [y, y, y], [y, y, z], [y, z, z], [z, z, z], 
  [x, x, x, x], [x, x, x, y], [x, x, x, z], [x, x, y, y], 
  [x, x, y, z], [x, x, z, z], [x, y, y, y], [x, y, y, z], 
  [x, y, z, z], [x, z, z, z], [y, y, y, y], [y, y, y, z], 
  [y, y, z, z], [y, z, z, z], [z, z, z, z]]

You can interpolate the matrix to turn it into a function over the reals rather than over the integers, which makes it easier to use with animate. In the example below, I made up a function to apply to a discrete dataset to simulate the two vectors of 2970 values that you're starting with.

restart:
N:= 2970:
stoimostyM:= sort(rtable(1..N, frandom(1..99, 1), datatype= hfloat)^+);
izmeniya:= (x-> sin(x^1.3)/sqrt(x))~(stoimostyM):
a:= <stoimostyM | izmeniya>:
f:= Interpolation:-LinearInterpolation(a):
plots:-animate(plot, [f, 1..X, numpoints= 500], X= 1..99, frames= 100);

Use polynomial arithmetic to construct all possible terms and then deconstruct the terms (without their coefficients) into lists by using op, `$`, and subsindets. Like this:

All:= (L::list(name), n::posint)->
    [subsindets(
        subsindets(([primpart]~@{op}@expand)(`+`(1,L[])^n - 1), `*`, op),
        `^`, `$`@op
    )[]]
:
All([x,y,z], 4);
[[x], [y], [z], [x, x], [x, y], [x, z], [y, y], [y, z], [z, z], 
  [x, x, x], [x, x, y], [x, x, z], [x, y, y], [x, y, z], 
  [x, z, z], [y, y, y], [y, y, z], [y, z, z], [z, z, z], 
  [x, x, x, x], [x, x, x, y], [x, x, x, z], [x, x, y, y], 
  [x, x, y, z], [x, x, z, z], [x, y, y, y], [x, y, y, z], 
  [x, y, z, z], [x, z, z, z], [y, y, y, y], [y, y, y, z], 
  [y, y, z, z], [y, z, z, z], [z, z, z, z]]

These are called delay differential equations. They can be solved numerically by Maple. See the help page ?dsolve,numeric,delay. You'll need to use the options delaymax and parameters.