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

ArrowOnTop:= proc(x::symbol, arrow::symbol:= `→`)
    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.

@Ahmed111 If we can assume that x, t, and lambda are real, then there'll be massive simplification. One way to do this is with evalc. Replace the command

simplify(subs(...), trig)

with

simplify(evalc(subs(...))) assuming (epsilon1, epsilon2, lambda1)::~complex;

Your simplify command begins

simplify((subs{

You need to change that to

simplify(subs({

The relation a <= b isn't exactly the same thing as the disjunction (a < b or a = b) because a <= b implies the existence of an order relation between a and b. In your example, the order relation wouldn't exist unless n were nonnegative. You can assume that this is so:

assume(n >= 0);

Then your context-menu-based test will return true.

If you use text-based commands, then an assuming clause can be used:

is(4/(9*n+7*sqrt(n)) <= 4/(9*n+7*sqrt(n))) assuming n >= 0;
                              true

The main difference between assume and assuming is that assumptions made with assume remain until they are removed, whereas those made with assuming are only in effect for a single command.

I posted my first Answer before you supplied an example procedure. It's not because your code is in a procedure that the Answer doesn't work; rather it's because of the combination of these two things:

1. The procedure is recursive, i.e., it calls itself;
2. The procedure produces its output with print statements rather than the normal "return value" mechanism.

There's nothing wrong with using recursive procedures; however, returning output with print (or any similar command) will usually lead to problems because you have no programmatic access to that output.

Here is your procedure:

Solveur := proc(b)
local a, e, i, k, L, B, H;
e := isolve(EQUATION2);
assign(e[1]);
a := {seq(x, k = 0 .. 10)};
L := [ ];
for i in a do if evalf(f(i)) > 1 then L := [op(L), f(i)];
end if; end do;
L; print(L);
unassign('x');
B := rhs(eval(e[1], _Z1 = 0));
if 0 < B then
Solveur(B);
end if;
end proc:

Presumably, EQUATION2 is some equation that you haven't supplied that contains x, b, and k. If that's not true, you'll need to provide a more-specific example. 

Change the procedure to this:

Solveur:= proc(b)
local k, a:= rhs(isolve(EQUATION2)[1]), B:= eval(a, indets(a, suffixed(_Z))=~ 0);
    [select(is@`>`, f~([{seq}(a, k= 0..10)[]]), 1)[], `if`(0<B, thisproc(B), [])[]]
end proc:

I can't test that without knowing your EQUATION2 and EQUATION1; however, I'll point out that neither the above nor your original could possibly produce 1 in their output because both check that f(i) > 1.

My procedure is meant to duplicate the results of your procedure except that the results are returned as a single list; however, if you wanted to generalize this for arbitrary EQUATION2 and EQUATION1, then there are several conditions that should be checked for, such as isolve not solving the equation or B not being numeric.

There are two things that I did that allowed for a massive simplification of the code:

1. I hard-coded t[j] = j/(M-1).
2. I generated the sequence of summation indices k[1], ..., k[s] by iterating over all size-s subsets of {0, ..., M-1} minus {j}. This works because the index notation that you supplied implies k[1] < k[2] < ... < k[s] (all inequalities strict).

My procedure below takes two arguments, t and M, and returns the entire list of polynomials l[j](t), 0 <= j < M. Please check that this gives the expected polynomials for some small values of M (because I may have over-simplified this).

Lpoly:= proc(t::algebraic, M::And(posint, Not(1)))
local
    beta:= proc(j::nonnegint, s::nonnegint)
    uses It= Iterator;
    local k, r, Mj:= Ms minus {j}, C:= (-1)^s/mul(j-~Mj);
        if s=0 then return C fi;        
        C*add(mul(Mj[k[r]+1], r= 1..s), k= It:-Combination(m,s))   
    end proc,
    m:= M-1, Ms:= {$0..m}, j, s
;   
    [seq](add(beta(j,s)*(m*t)^(m-s), s= 0..m), j= 0..m)
end proc
: