The command is plots:-display. See the help page ?plots,display.

Here's some code for it:

restart:
f:= (x,y)-> 100*(y-x^2)^2 + (1-x)^2: #objective
fx:= D[1](f);  fy:= D[2](f); #its gradient
(initx, inity):= (-1, 1):

PathMin:= proc(a,b)
local 
    t, #unassigned symbolic

    #x- and y- coordinate functions along steepest-descent line:
    X:= unapply(a-fx(a,b)*t, t),
    Y:= unapply(b-fy(a,b)*t, t),

    #objective function restricted to steepest-descent line: 
    P:= unapply(f((X,Y)(t)), t),

    R:= [fsolve](D(P)(t), t= 0..infinity),
    m:= min(P~(R))
;
    if R::{[], specfunc(fsolve)} then
        print("No critical points along line");
        return
    fi;
    #new values for a and b:
    (X,Y)(select(r-> P(r)=m, R)[1])
end proc
:
(a1,b1):= (initx,inity):
(a1,b1):= PathMin(a1,b1);
                    a1, b1 := 1.000000000, 1

So, we got to the global minimum in one step. I think that's just due to having a super-lucky choice of initial point. 

We can generalize from 2*Pi/14 to 2*Pi/n, use symbolic summation, and still get the simplification. Like this:

restart:
w:= 2*Pi/n:
v:= <1, (sin,cos)(w*t)>:
simplify(sum(v.v^%T, t= k..k+n-1));
             

Note that the LinearAlgebra package and the commands with, Vector, and Transpose are not needed. The notation (...)^%T is a transpose operator; it can also be (...)^+ in 1D (plaintext) input. The <...> is a column-vector constructor. The (f,g)(x) is equivalent to f(x), g(x).

If your report is accurate, then that's a bug, but that bug has long since been fixed. I see that you're using Maple 18. In Maple 2019, I get

Logic:-Dual(a &implies b);
               `&not`(b &implies a)

And that is equivalent to `&not`(a) &and b.
 

Once you have the module pds, you can do this to be reminded of the module-specific commands available (which are called the module's exports):

exports(pds);
             plot, plot3d, animate, value, settings

Some examples of using those:
pds:-plot3d(u(x,t), x= 0..1, t= 0..1);

pds:-plot(theta(x, t), x= 0..1, t= 0.1);

pds:-value(u(x,t))(0.2, 0.1);
      [x = 0.2, t = 0.1, u(x, t) = 0.0597342089921277242]

Here's a procedure for it. No package is needed:

ArrowOnTop:= proc(x::symbol, arrow::symbol:= `&rightarrow;`)
    nprintf(`#mover(mn(%a), mo(%a))`, cat("", x), cat("", arrow))
end proc
:
ArrowOnTop(x);
             

To get the x in italic, change mn to mi.

In most or all of the cases that you're interested in, a formula is known for arbitrary terms of the power series, and that formula can be presented as an inert summation like this:

'tan@@(-1)'(x) = convert(arctan(x), FormalPowerSeries);

FormalPowerSeries can be abbreviated FPS.

So, this doesn't use the series command nor the series data structure at all (not even in the background) and thus completely gets around the issue of O(...) or `...`.

Here I present six one-line methods to do it. First, I define some variables for the sake of generalizing your example:

d:= 4:               #degree
M:= 3:               #max coefficient + 1
R:= [M$d-1, M-1];    #radices
j:= [d-k $ k= 1..d]; #exponents
A:= a||~j;           #symbolic coefficients
                       R := [3, 3, 3, 2]
                       j := [3, 2, 1, 0]
                     A := [a3, a2, a1, a0]

Note that the product of the radices is 3*3*3*2 = 54.

Each of the six methods below generates a list (of 54 elements) of lists (each of 4 elements) of coefficients. In a later step, they'll be made into the polynomials, and 1+x^4 will be added to each. Two of the methods use an embedded for-loop; these will only work using 1D input (aka plaintext or Maple Input). Three use seq for the loops, and they'll work in any input mode. One uses elementwise operation (~) instead of a loop.

The first two methods use an Iterator called MixedRadixTuples. This is very similar to the CartesianProduct shown in VV's Answer. MixedRadixTuples is specifically a Cartesian product of sets of consecutive integers, each set starting at 0. The only difference in the two methods below is that the second uses seq instead of for.

C1:= [for C in Iterator:-MixedRadixTuples(R) do [seq](C) od];
C2:= [seq]([seq](C), C= Iterator:-MixedRadixTuples(R));

The next three methods generate the lists of coefficients from the base-3 representations of each integer from 0 to 53.

C3:= [for m to mul(R) do convert(M^d+m-1, base, M)[..-2] od];
C4:= [seq](convert(M^d+m-1, base, M)[..-2], m= 1..mul(R));
C5:= convert~(M^d+~[$0..mul(R)-1], base, M)[.., ..-2];

The final method generates a nested seq command and then executes it.

C6:= [eval](subs(S= seq, foldl(S, A, seq(A=~ `..`~(0, R-~1)))));

Now I show that each of the above lists is the same by showing that the set formed from all six has only one element.

nops({C||(1..6)});
                               1

The final list of polynomials can be formed via a matrix multiplication, the coefficients as a 54x4 matrix times a column vector of the powers of x:

po:= [seq](rtable(C1).x^~<j> +~ (1+x^d));

A variation of that last command is

po:= map(inner, C1, x^~j) +~ (1+x^d);

The undocumented long-standing built-in command inner does the dot (or inner) product of two lists as if they were vectors.

The following command returns a table whose indices are the output of galois and whose entries are the corresponding subsets of the input list s:

GG:= ListTools:-Classify(galois, select(irreduc, s)):

If you then want to count those subsets, do

GGn:= nops~(GG);

If you just want one of the five given representations of the groups, then replace galois in the first command with p-> galois(p)[k], where 1 <= k <= 5.

Or the desired information can be extracted from table GG like this:

rtable([curry(op,[2,-1])@[lhs], nops@rhs]~([indices](GG, pairs)));

Here's how to make the loop-based methods (either for or seq) efficient. The key thing is to not use a list for S while the computations are being done; instead, make S a sparse table. A table is one of the fundamental container objects in Maple (along with lists, sets, and rtables). If a table is declared sparse, then unassigned entries evaluate to 0. The table can be efficiently converted to a list when the computations are done. Here's your code (in the box below). Code in green can remain unchanged; that in red must be changed.

roll:=rand(1..9):
N := sort([seq(seq(10*(2+jdx)+roll(),jdx=1..5),idx=1..10)]):
V := 10*seq(roll()+10*trunc(N[idx]/10),idx=1..nops(N)):

## generate index for sum. This groups the values for sum
sumidx := [seq(floor(N[idx]/10)-2,idx=1..nops(N))];

## create and zero sum variable
S := 0 *~convert(convert(trunc~(N/10),set),list); #Corrected in note 1 below

## calc sum and display it - This works
for idx to nops(sumidx) do
    S[ sumidx[idx] ] += V[idx]:  
end do:
S; #See note 2.

## error because sumidx[idx] is not evaluated ?????
S := 0 *~convert(convert(trunc~(N/10),set),list): #See note 1. 
seq('S'[ sumidx[idx] ] += V[idx], idx=1..nops(N)); #See note 3. 
S; #See note 2.

And here are the notes with my corrections:

1. In both cases, the line should be replaced with
S:= table(sparse):
There's no need to create a table with any specific number of elements; indeed, there's not even a formal mechanism for doing that if you wanted to.

2. In both cases, the line should be replaced with
S:= [entries](S, indexorder, nolist);
This is the efficient conversion of the completed table into a list. A commonly used more-intuitive variation that's not quite as efficient but still fine is
S:= convert(S, list);

3. The correction of this line has already been the subject of several Answers and Replies, and you even had a correction in your Question. My preference is 
seq((S[sumidx[idx]]+= V[idx]), idx= 1..nops(N)):
although the differences between that and your f or acer's `+=` are minor.

So, you can see that converting accumulator variables from lists to tables is quite simple. If the lists are long, the time savings will be huge.

Here's another easy adjustment to correct your original. The reason that your seq doesn't work is because embedded assignments, including embedded updating assignments such as +=, always need to be in parentheses (even when they can't possibly disambiguate anything). So, you just need to correct your line

seq('S'[ sumidx[idx] ] += V[idx], idx=1..nops(N));

to

seq( ( S[sumidx[idx]] += V[idx] ), idx=1..nops(N));

And you need to remove the unevaluation quotes. Actually, I don't think that they ever did anything useful in these examples.

All four of the following methods are atrociously inefficient (and all for the same reason): this Answer, your for-loop method, your seq using procedure f, and acer's seq using `+=` (but I know that he only did that because he wanted to show you a minimal correction to what you already had, so I'm not criticizing). The reason they're inefficient is that they make assignments to list elements. Every iteration recreates the entire list S even though only one element is being changed. I'm writing improvements to these loop-based methods that avoid this problem, which I'll post soon.

[Edit: While adapting the OP's code to my methods, I initially missed his 10 immediately preceding trunc in the construction of V. I just added it to the code below, before the iquo, which I replaced trunc with.]

[Edit 2: I also mixed up the roles of rows and columns in the matrix, and that's now corrected below.]

It can be done much more simply than your original by making V a 5x10 matrix and then adding its rows. Like this:

restart;
roll:= (n::posint)-> RandomTools:-Generate(list(posint(range= 9), n)):
nN:= (Jx:= 5)*(Ix:= 10):  
N:= sort(10*['seq'(3..2+Jx)$Ix] + roll(nN));
V:= Matrix((Jx,Ix), 10*(roll(nN) + 10*iquo~(N,10)));
S:= [seq](ArrayTools:-AddAlongDimension(V,2));

The final output is now 
        S := [3420, 4570, 5480, 6620, 7320]
matching yours and acer's S.

This technique requires no indices whatsoever.

It can be done with a single plot command and directly nested seqs, by which I mean one seq command immediately inside another with no command between them. Like this:

restart;
a:= (k,p)-> (100-p)*exp(-k/2):
pvals:= [seq](0..200, 50):
colors:= [blue, red, green, magenta, cyan]:
opts:= 
    style= point, symbol= solidcircle, symbolsize= 15,
    color= colors, legend= (p=~ pvals)
:
plot([seq]([seq]([k, a(k,p)], k= 0..10), p= pvals), opts);

So all these are not needed: the plots package and the with, display, and pointplot commands.
 

ModuleIterators

Also, it is possible to have a single seq with multiple indices, although for this particular plotting example it's not as useful as the nested seq above. Here's an example:

opts:= style= point, symbol= solidcircle, symbolsize= 15:
plot(
    [seq]([j[1], a(j[1],j[2])], j= Iterator:-CartesianProduct([$0..10], pvals)),
    opts
);

This type of index is called a ModuleIterator, and they can be used as indices for seq, add, mul, etc., commands and for-loops. You can write your own ModuleIterators, and numerous pre-written ones are provided in the Iterator package.

I am very skeptical of your time measurements because they aren't strictly increasing. You may not be properly accounting for garbage collection, for example. For each number of digits, take an average for several cases with different input.

You should only use strictly increasing model functions, so no polynomials with more than one term.

The best fits that I got for your data were with the model function T(n) = a*n^b*c^n, which can be linearized as ln(T(n)) = ln(a) + b*ln(n) + ln(c)*n.

dofit := proc(N::list(numeric), T::list(numeric), n::name)
local O:= exp(Statistics:-LinearFit([1,ln(n),n], N, ln~(T), n));
    print(plot([`[]`~(N,T), O], n= .9*N[1]..1.1*N[-1], style= [point,line]));
    simplify(O) assuming n::positive
end proc
:

dofit(N[..-2], Tcpu, n);

I didn't take the time to find the exact cause of your error. Rather, I just quickly changed everything that was even remotely suspicious to me, and it worked. I got

restart;
N := 2:
var:= seq(alpha[n], n = 1 .. N):
v:= add(alpha[n]*phi[n](x), n = 1 .. N): #Use add, not sum
R:= diff(v, x $ 4) - p/EI:
eqs := seq(int(R*phi[n](x), x = 0 .. L) = 0, n = 1 .. N);
phi:= x-> sin((2*op(procname) - 1)*Pi*x/L); #op(procname) refers to index n.      
solve({eqs}, {var});

The op(procname) in the definition of phi represents its index, n. This crude mechanism is the only way (I know of) to use parameters passed as indices rather than as arguments.