@Carl Love I agree you. Yet I cannot find any detailed explanation that "Maple implementation is very inefficient" in the official documentation. (In other words, what I need is not simply a known fact or practical experience.) 

Actually, I read some viewpoints like: GOTOs can be a useful language feature improving program speed, size and code clarity, and instructions like break and continue are just the old goto in sheep's clothing. Is there any difference between Maple's goto and FORTRAN/C/...'s goto?

Another statement (from the Wikipedia article GOTO) is: Functional programming languages such as Scheme generally do not have goto, instead using continuations. But what about Maple?

@C_R You can change the third line to 

isolate(DETools:-firint(ode*DEtools:-intfactor(ode, phi(s)), phi(s)), _C1); # instead of c__1, but I don't know why the original one works before

As for the second error message, the ?? should be (1) (which is an equation label).

@ Sorry. What does "previous and more robust" mean? If I understand right, do you mean that the current version is less robust instead? 

Besides, it seems that neither the built-in Optimization package nor the additional GlobalOptimization package is well maintained; I cannot find any relevant news in most recent "Release Notes". (Actually, in some cases, MatLab's Optimization Toolbox and Global Optimization Toolbox are more powerful, but they are not easy-to-use for advanced users that are just familiar with Maple syntax.)

@Carl Love Amazing! Your code works well. Many thanks!
In addition, in that link, the (Mathematica) example uses 

but it appears that there is not a similar built-in functionality (Am I wrong?) in Maple. So, to verify the result, I have to copy them to Maple (and then execute `convert/FromMma`), which is more or less inconvenient. Since I never want to make an ugly procedure again, can this be implemented in a single line of input elegantly?
Unfortunately, the following implementation is quite inelegant, and does not work in general cases: 

(*  list all possible trees formed by certain binary combiners inefficiently, where leaves are distributed keeping the order at the deepest level (if possible, as the result below goes wrong) *)
(z -> ListTools:-FlattenOnce(((y -> ListTools:-MakeUnique(ListTools:-FlattenOnce(`~`[x -> ['[x[() .. k - 1][], F(x[k], x[k + 1]), x[k + 2 .. ()][]]' $ ('k' = 1 .. numelems(x) - 1)]](y))))@@(numelems(z[]) - 1))(z)))([[a, b, c, d, e]]);
 = 
    [F(F(F(F(a, b), c), d), e), F(F(F(a, b), c), F(d, e)), 

      F(F(F(a, b), F(c, d)), e), F(F(a, b), F(F(c, d), e)), 

      F(F(a, b), F(c, F(d, e))), F(F(F(a, F(b, c)), d), e), 

      F(F(a, F(b, c)), F(d, e)), F(F(a, F(F(b, c), d)), e), 

      F(a, F(F(F(b, c), d), e)), F(a, F(F(b, c), F(d, e))), 

      F(F(a, F(b, F(c, d))), e), F(a, F(F(b, F(c, d)), e)), 

      F(a, F(b, F(F(c, d), e))), F(a, F(b, F(c, F(d, e))))]

Any ideas?

@nm Actually, Mma's LeafCount has an alternate (open-source) official version (though this may not give the total number of indivisible subexpressions): The automatic complexity function.
Besides, Maple's dismantle should be used when the "size" is measured (since the `rtable`, especially the `Matrix`, is not equivalent to nested lists). However, it's hard to sort a list of expressions by Depth and then LeafCount likewise in Maple.

@C_R I find a way to "simplify" these, though it is not guaranteed in general. Fortunately, it works here: 

alias(seq(alpha[i] = RootOf(_Z^3 + _Z - 1, index = i), i = 1 .. 3)):
f__5 := collect(evala(collect(f__0, {y, x}, 'distributed', n -> convert(RootOf(PolynomialTools[MinimalPolynomial](evalf(n), _X, 3), _X, evalf(n)), RootOf))), x, evala);
 = 
[ 3   1 /        2                 \          
[x  + - \alpha[1]  + 2 alpha[1] + 2/ (y + 2) x
[     3                                       

     1 /          2               \ / 3                     \   3
   - - \2 alpha[1]  + alpha[1] + 4/ \y  - 2 alpha[1] + y + 1/, x 
     3                                                           

     1 /        2                 \          
   + - \alpha[3]  + 2 alpha[3] + 2/ (y + 2) x
     3                                       

     1 /          2               \ / 3                     \   3
   - - \2 alpha[3]  + alpha[3] + 4/ \y  - 2 alpha[3] + y + 1/, x 
     3                                                           

     1 /        2                 \          
   + - \alpha[2]  + 2 alpha[2] + 2/ (y + 2) x
     3                                       

     1 /          2               \ / 3                     \]
   - - \2 alpha[2]  + alpha[2] + 4/ \y  - 2 alpha[2] + y + 1/]
     3                                                       ]


verify(p, -3*`?()`(`*`, f__5), equal);
 = 
                              true

Compare:  

> length[shortname](f__5);
                                                          771

> MmaTranslator[Mma][LeafCount](f__5);
                                                          247

@acer You are right. However, if one can obtain an equivalent but still readable (exact) absolute factorization of p (I mean, in another way), this is also welcome. 
Edit. I find a Magma program to do so: Programmes Magma, Maple.

@Carl Love The reason is: 

It appears that the former one does more things (e.g., prevents premature evaluation, as described in this recent question) internally, hence in my view, `$` should be a low-level level. 

@Derek Wright Thanks for your concise comparisons. Unfortunately, this link appears dead. Are such comparisons still available? By the way, it seems that Maple's numerical solvers have certain limitations on solving differential equations: A Comparison Between Differential Equation Solver Suites In MatLab, R, Julia, Python, C, Mathematica, Maple, and Fortran - Stochastic Lifestyle. Will this degrade simulation accuracy and even cause major errors in some model?

@Harry Garst In my opinion, seq should also work here, but I found some issues before, so I have to use the lower-level command `$` at last. By the way, how about the so-called Chió's condensation (compared to Dodgson condensation)?

@ Well, I'm not sure if Maple includes (optimized) implementations of those algorithms from recent research which routinely outperform the classic and standard C/Fortran methods. Any benchmarks (liike Benchmarks for SciML and Equation Solvers)? 

I don't know what "tolerance = tolerance" means. If you use "res := eval(<vars>, fsolve(eval({e||(1..8)}, [y1[n]=iny1, y2[n]=iny2]), {vars})):" instead, your code will work well.

Not an answer. But if you use "[x]" (instead of "x" or "{x}"), the output does agree with each other: 

with(RealDomain):
sol1 := solve(x - 1 > 0, [x]);
sol2 := solve(1/x^2 > 0, [x]);
sol3 := solve(1/(x^2 - 1) > 0, [x]);
sol4 := solve(x^2 + 1 > 0, [x]);
sol5 := solve(x^2 + 1 > 0, x); = 
                       sol1 := [[1 < x]]

                       sol2 := [[x <> 0]]

                  sol3 := [[x < -1], [1 < x]]

                       sol4 := [[x = x]]

                           sol5 := []

@acer Thank you. Though this might produce useful solutions (especially in physics and engineering), I think that for the original problem, this does not guarantee the completeness of results.