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
:    

I think that this is the easiest way (and it's also very efficient), but it only works since Maple 2018 and only using 1D input (plaintext input). Let's say that you have a loop that's producing output with each entry on its own line, such as

for k to 9 do
    ithprime(k)
od;

Then change that to

L:= [
    for k to 9 do
        ithprime(k)
    od #no semicolon here!
];

That works for any do loop; it needn't be a for-do loop.

The line-breaks and other spacing in my code are just for emphasis. You could just as well use

L:= [for k to 9 do ithprime(k) od];

If the original loop delivers its output via a print command, then you'll need to remove the print.

Only the values produced on the last "line" of the loop (i.e., right before the od) will be included in the list.

The analytic solution for x(t) in the ODE system in your original Question is the constant function x(t)=1. (This is obvious, but let me know if you need further explanation of why that is.) In your followup, you plotted the analytic solution functions x(t) and y(t) for a similar---but not identical---ODE system. The crucial difference is that x(t) is not constant there because you changed its initial conditions.

The output of odeplot(dsolve(..., numeric)) is not a parametric plot of y(t) vs. x(t); it is a plot of x(t) vs. t. In that plot, there is some miniscule variation of x(t) due to the numeric approximation process used by dsolve. This variation is entirely in the 16th decimal place and beyond, which makes it many orders of magnitude smaller than the stated error limits for numeric dsolve. Since x(t) is otherwise constant, this variation in the only variation in the vertical coordinate, so it alone determines the scaling of the vertical axis. It is unreliable and inefficient to try to control this with numpoints (won't work) or Digits (might or might not work; inefficient). If for some reason you actually want a plot whose only content is a constant function, use the view option to set the scaling (for example, view= [0..10, 0.9..1.1]).

Here's how to tell odeplot the functions that you want to plot:

restart;
g:= 9.81;
sol:= dsolve(
    {
        diff(x(t), t, t) = 0, diff(y(t), t, t) = -g + y(t)^2, 
        x(0) = 1, y(0) = 0, D(x)(0) = 0, D(y)(0) = 0
    }, 
    numeric
):
#Plot x(t) and y(t) vs. t on the same axes:
plots:-odeplot(sol, [[t,x(t)], [t,y(t)]], t= 0..10); 

It doesn't make sense to do a parametric plot in this particular case, one coordinate being constant. But otherwise the function specification could be [x(t), y(t)].

If you make the function specifcation [t, x(t) - 1], you'll see that all the variation in x(t) is in the 16th decimal place and beyond.

That should be 

f:= x-> piecewise(...)

just like you had in your earlier Questions. (I wonder why you didn't notice that?)

That's quite a weird error message though! Don't bother trying to understand it at this point; I think that it may be intended primarily for developers.

The command lprint(eval(r)) could be used as an aid for checking whether you've entered r as intended. All implicit multiplications will be shown explictly as *. There are two spots where I think that you intended multiplication, but it doesn't appear, which I've highlighted below with ^^.

lprint(eval(r));
(phi, B, theta) -> sqrt(-B^2+(a^2+B^2)*cos(theta)^2*f)(sqrt((-B^2+(a^2+B^2)*cos
                                                     ^^
(theta)^2)/(a^2+B^2)/cos(theta)^2)(phi-sqrt(a^2+B^2)*cos(theta)/sqrt(-B^2+(a^2+
                                 ^^
B^2)*cos(theta)^2)*arcsec(a/sqrt(-B^2+(a^2+B^2)*cos(theta)^2))))

2D Input usually requires a space between operands in order for it to be interpreted as multiplication. This is especially important when the right operand is in parentheses because the syntax checker usually won't catch the error in that case (because such expressions have a legitimate interpretation other than multiplication).

Also, the lprint shows that f is the last operand "under" the first square root. Did you intend for f to be outside that square root instead?

If I make these changes, I get this plot:

Is the returned p a local variable, perhaps local to the procedure that returned it? That would be the most likely explanation. 

Otherwise, we need to see the code that returns p.