@vv Thank you. Unfortunately your code does not work on my computer for some reason:

Only with the option  compile=false  I get the result, but the time is 10 times longer.

When I wrote about my method for 100 seconds, I used another technique (uncompiled code). I'll submit this code later.

@tsunamiBTP  I did this and did not find any fundamental discrepancies. Some differences are due to the fact that in  M_new  I am greatly rounding the results (evalf[3]). For comparison, I calculated the hundredth line of my and your matrix. The graphs are indistinguishable. See untitled5_new2.mw

@John Fredsted  About this construction see Matrix and Vector Construction Shortcuts help page. By the same method, we can remove from a matrix any row or several rows or columns:

A:=LinearAlgebra:-RandomMatrix(6, generator=0..9);
A1:=<A[..3], A[5..]>;  # 4th row deleted of A
A2:=<A[..,..4] | A[..,6]>;  # 5th column deleted of A

I think this method is faster than calling  LinearAlgebra:-DeleteRow  or  LinearAlgebra:-DeleteColumn .

@tsunamiBTP  See corrected file  untitled5_new1.mw 

@tsunamiBTP  You wrote "This seems to work like a nested FOR loop". No, these are the different things. Here is a simple example of use nested seq. Let's say you want to create a matrix 10x10 that has symbols  x1, x2, ... ,x10  in the first line, the second line contains the same symbols, which squared, then cubed and so on:

Matrix([seq([seq((x||j)^i, j=1..10)], i=1..10)]);

Another example. It is very convenient and effective to use  seq  command  to obtain the sequence or list of objects with specified properties. In the example, we generate the sequence of Pythagorean triplets with  hypotenuses <=100:

seq(seq(seq(`if`(a^2+b^2=c^2,[a,b,c],NULL), b=a+1..c-1), a=1..c-2), c=3..100);

[3, 4, 5], [6, 8, 10], [5, 12, 13], [9, 12, 15], [8, 15, 17], [12, 16, 20], [7, 24, 25], [15, 20, 25], [10, 24, 26], [20, 21, 29], [18, 24, 30], [16, 30, 34], [21, 28, 35], [12, 35, 37], [15, 36, 39], [24, 32, 40], [9, 40, 41], [27, 36, 45], [14, 48, 50], [30, 40, 50], [24, 45, 51], [20, 48, 52], [28, 45, 53], [33, 44, 55], [40, 42, 58], [36, 48, 60], [11, 60, 61], [16, 63, 65], [25, 60, 65], [33, 56, 65], [39, 52, 65], [32, 60, 68], [42, 56, 70], [48, 55, 73], [24, 70, 74], [21, 72, 75], [45, 60, 75], [30, 72, 78], [48, 64, 80], [18, 80, 82], [13, 84, 85], [36, 77, 85], [40, 75, 85], [51, 68, 85], [60, 63, 87], [39, 80, 89], [54, 72, 90], [35, 84, 91], [57, 76, 95], [65, 72, 97], [28, 96, 100], [60, 80, 100]

@Markiyan Hirnyk  Here is a check (I just added a few commands to your document):

Download Check.mw

@mehdibaghaee Sorry, I did not notice the second condition.
Should be:

seq(`if`(V[i]>0 and V[i]<>infinity, [i, V[i]], NULL), i=1..10);

@Preben Alsholm

For example, if we have a function of one variable, given explicitly  f:=x->f(x)  and we want to  calculate the derivative value at a point  x=x0 , then I think that  D(f)(x0)  is better than  eval(diff(f(x),x), x=x0)  or  unapply(diff(f(x),x), x)(x0) .

@Preben Alsholm 

D operator can be used in similar situations, only the function must be specified explicitly:

restart;
f:=(x,y)->(x^3+y^3)^(1/3);
r:=(s,t)->(s*sin(t), s+56);
D[1](unapply((f@r)(s,t), s,t));
 

@tsunamiBTP

1. Since you mentioned  plottools:-getdata  command, I understood your request as a desire to get data from the plot.

2. If you want to get data from the second your graph, write:

data := op([2, 1], P);  # data as a matrix
convert(data, listlist);   # data as a list of lists

@Carl Love  For greater clarity, you can add one row to the final matrix:

<< x(`in degrees`) | sin(x) | Error> ;  < d | S | Err >>;

@tomleslie   I understood this condition  0<= beta, delta <= 1, that both roots should be searched in the range 0..1. If we follow your understanding, then here are 4 plots (I used to believe my eyes). The first two graphs clearly show the roots that are not in the list of DirectSearch. The last 2 plots show that the last 2 solutions of DirectSearch are wrong.

plots:-implicitplot([focdeltapioptS2Tbeta_eg, focbetapioptS2Tbeta_eg], beta=0..2, delta=-20..1, color=[red,blue], thickness=2,  gridrefine=5, axes=normal);

plots:-implicitplot([focdeltapioptS2Tbeta_eg, focbetapioptS2Tbeta_eg], beta=0.99..1.01, delta=-50..1, color=[red,blue], thickness=2,  gridrefine=5, axes=normal);

plots:-implicitplot([focdeltapioptS2Tbeta_eg, focbetapioptS2Tbeta_eg], beta=10.8..11, delta=-2.2..-2, color=[red,blue], thickness=2,  gridrefine=5, axes=normal);

plots:-implicitplot([focdeltapioptS2Tbeta_eg, focbetapioptS2Tbeta_eg], beta=14.7..14.9, delta=0.9..1.1, color=[red,blue], thickness=2,  gridrefine=5, axes=normal);
 

@Wtolrud  You can use a simple procedure named  Sort  that sorts any polynomial in ascending order of degrees. Formal parameters of the procedure: P is a polynomial (possibly with symbolic coefficients), var is a polynomial variable (by default this is x):

Sort:=proc(P::polynom, var::name:=x)
sort(P, var, ascending);
end proc:

Examples of use:

Sort(a*x^2+b*x+c);
Sort(-3*t^3+t+2*t^2+5, t);

                                   c+b*x+a*x^2
                                 5+t+2*t^2-3*t^3
 

@omkardpd  Place your specific system here in a copyable form as a text (not as a picture).

You forgot the sign of multiplication. Should be:

5*x^3+7*x^2+2*x^3