@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

@Ramakrishnan  I took frames = 121 so that one of the animation frames matched the value of the parameter a = -3, at which a qualitative change in the structure of the solution occurs. Therefore, the step for the parameter  a  can be taken from several variants: 0.5, 0.2, 0.1, 0.05, ... . I chose the latter option and then the number of frames will be  6 / 0.05 + 1 = 121

@Ramakrishnan  Substantially reduce the number of frames, for example to 60, or write so:

restart;
A := plots:-animate(plot, [x^3+a*x+2, x = -4 .. 4, -15 .. 15, color = blue, thickness = 2], a = -6 .. 0, frames = 120): 
A;

Also note that I slightly simplified the penultimate line of my code above (for C).
I advise you to think about why I took in the code  frames=121  instead of frames=120 .

@vv  Very interesting how did you get these results? My way for some reason fails for large n (even option remember does not help), and Markiyan's way takes a long time:

P(10^5);
   Error, (in P) too many levels of recursion

restart;
ts:=time():
x(1) := 1: N := 100000; for n to N do x(n+1) :=
sum(convert(x(n), base, 10)[j]*10^(nops(convert(x(n), base, 10))-j), 
j = 1 .. nops(convert(x(n), base, 10)))+n end do:

x(10^5);
time()-ts;

                          N := 100000
                           9966045232
                             84.281
 

@Markiyan Hirnyk  I often visit this site  http://www.mathforum.ru/forum/list/1/

Your expression can be greatly simplified if you reduce the fraction by  sqrt(s+thetac)  because

         A1*Dc*alpha1*s^2+A1*Dc*alpha1*s*thetac=A1*Dc*alpha1*s*(s+thetac)

@MrMarc  If you need a single solution, then you need to impose an additional condition. For example, the condition  x=y  leads to the unique solution:

solve({x+y = 373320, z = (x+y) / 0.44 - y -  y* (1 - 0.99), x=y});

         {x = 186660., y = 186660., z = 659927.9456}
 

 rlopez Many thanks for the detailed analysis of the code in  IntOverDomain  procedure and the great work you have done for this.