@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.

The spectral radius is also the maximal singular value^[*1]. Regardless of whether you use command Eigenvalues or SingularValues from LinearAlgebra, numeric results can be quickly obtained by converting the Matrix to datatype= hfloat.

restart:
SpectralRadius:= proc(G::Graph)
uses LA= LinearAlgebra, GT= GraphTheory;
    (max @ LA:-SingularValues)(
        rtable(GT:-AdjacencyMatrix(G), 'datatype'= 'hfloat')
    )
end proc
:
GT:= GraphTheory: SG:= GT:-SpecialGraphs:
M22:= SG:-M22Graph(); #This graph is 16-regular.
      Graph 1: an undirected unweighted graph with 77 vertices 
      and 616 edge(s)

CodeTools:-Usage(SpectralRadius(M22));
memory used=241.62KiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=6.00ms, gc time=0ns

                        16.0000000000000

[*1] The singular values of a real or complex matrix A (not necessarily square) are the square roots of the eigenvalues of A . HermitianTranspose(A). These are necessarily nonnegative reals. In the square floating-point case, these can be computed more efficiently and more stably than the eigenvalues of A itself.

The _Y could be replaced by any otherwise unused variable. It's what logicians and computer scientists call a bound variable, and here it's bound to the DESol expression. Mathematicians sometimes call it a dummy variable, but I disdain that term. So, the important thing is what DESol means rather than what _Y(x) means.

So, Maple can't solve your ODE, but it can simplify it a tiny bit. It's saying, essentially, "Here's another, perhaps simpler, ODE (the 1st argument of the DESol equated to 0) whose solution is the same as that of your original ODE." Now, you may think that this new ODE is such a trivial change from your original that it's hardly worth stating, and if all DESols returned by dsolve were that trivial, then I'd agree. But sometimes they're not so trivial.

You could try exporting the worksheet as PDF and then copy-and-paste from the PDF. From the worksheet's File menu, select Export As, then select PDF from the "Files of type" pull-down.

A variation of Tom's Answer:

select(j-> j[[4,7]] = [[C,Z], [B,Y]], Perms)

One way:

solve(2*a +~ (b+2)*~[0,1] =~ [1,1]);

Another (this one doesn't use elementwise operations):

x:= t-> 2*a+b*t+2*t:  solve({x(0)=1, x(1)=1});

Another:

solve((t-> 2*a+b*t+2*t)~([0,1]) =~ [1,1]);