@Carl Love  I meant that when calculating a single expression  Expr  (if we need  n  digits) it is not practical to reset the value of  Digits . We can just write 
evalf[n](Expr)  or  evalf(Expr, n)

@asa12  I noticed that you are using Maple 12, which does not have  numelems  command (it was introduced in Maple 15) and also there are no combinat[nextperm], combinat[prevperm], combinat[firstperm] and combinat[lastperm] commands (they were introduced in Maple 16). Instead of  numelems  command, for lists you can use  nops  command, and for getting the needed permutations you can use  Iterator  package (Joe Riel already wrote about it above), which you can freely download from the Maple application center from this page 

x:=infinity*I;
x^2;
exp(%);
 

@Carl Love  Thanks. I have never used  Grading  package before. I just found a suitable example in the help on  examples,Quiz (Quiz)  and slightly changed it.

@minhthien2016 

restart;
k := 0:
for a from -5 to 15 do
for b from -5 to 15 do
for c from -5 to 15 do
if igcd(b, c) = 1 and nops({a, b, c}) = 3 then X := [solve(x = a+sqrt(b*x+c))];
if nops(X)=2 and type(X[1], integer) and type(X[2], integer) then k := k+1; L[k] := [a, b, c, X[]] end if end if end do end do end do;
L := convert(L, list):
k;
L;

@minhthien2016  I found an error in my code above . Here is the correct code:

restart;

k:=0:
for a to 10 do
for b to 10 do
for c to 10 do
for d to 10 do
for e to 10 do
for f to 10 do

if a*e-b*d > 0 and a^2*f^2-a*b*e*f-2*a*c*d*f+a*c*e^2+b^2*d*f-b*c*d*e+c^2*d^2 < 0 and igcd(a, b, c) = 1 and igcd(d, e, f) = 1 and igcd(a, b, c, d, e, f) = 1 then X:=[solve(x = (a*x^2+b*x+c)/(d*x^2+e*x+f))] ; 
if type(X[1],integer) and type(X[2],integer) and type(X[3],integer)  then k:=k+1;  L[k]:= [a, b, c, d, e, f, X[]] fi; fi;

od: od: od: od: od: od:

L:=convert(L, list):
 
k; L;
      

@Carl Love  Yes, of course you are right. Should be

plots:-implicitplot([x^2*y^2-2*x*y^3+y^4-y^3+x^2-1=0, x=y+sqrt(y), x=y-sqrt(y), x=y+sqrt(y-1), x=y-sqrt(y-1)], x=0..10, y=0..13, color=[red, blue$2, green$2], gridrefine=4, size=[1000,1000]);

Thank you.
 

@Carl Love  For one branch it is closer, for another - farther. Take a look:

A:=plots:-implicitplot(x^2*y^2-2*x*y^3+y^4-y^3+x^2-1=0, x=0..10, y=0..13, color=red, gridrefine=4):
B:=plot([y+sqrt(y), y-sqrt(y)], y=0..10, color=blue):
C:=plot([y+sqrt(y-1), y-sqrt(y-1)], y=0..10, color=green):
plots:-display(A, B, C);

I used the following:

solve(x^2*y^2-2*x*y^3+y^4-y^3+x^2-1=0, x);
map(asympt, [%], y, 0);
 

@basha 666  I showed in principle for your average sample  S=0  how it is possible to plot contour lines with the filling of individual regions by the specific colors. In my example, the program selects contour lines automatically (in the code, only their number is indicated). If you want to plot contours for some specific values  C  (here  psi(eta, x)=C ), you must specify these values  C .

@tomleslie In the text part of the help on plots:-arrow, the option  color  is not really mentioned, but in the examples it is used. Therefore, it is natural to think that it was simply forgotten to mention in the text part of the help. On the other hand, this option is also absent in the text part of  plots:-display  command.
If we follow your recommendation, then in order to plot even one arrow with the needed color, we should call plots:-display  command?...

@ssllys  I just noticed that  af_eq  is quadratic with respect to sigma. Solving this quadratic equation, we get 2 branches  curve1  and  curve2  (right and left). Here is an example of working with these branches  (I plot these branches in different colors):

af_eq:=(3*alpha[2](t[1], t[2])*A^3*(1/8)-A*sigma)^2+(1/4)*alpha[1](t[1], t[2])^2*A^2 = (1/4)*f0^2;
f0 :=0.1; alpha[1]:=0.0087;alpha[2]:=2.5871;
af_eq;
curve1,curve2:=solve(af_eq, sigma);
minimize(curve1, A=0.1..0.5, location);
A0:=rhs(%[2,1,1,1]);
with(plots):
P1:=implicitplot([sigma=curve2,sigma=curve1],sigma = -1.2 .. 2, A = 0 .. A0,  numpoints = 20000, axes = box, axesfont=[SYMBOL, 14],labels = [sigma, A], labelfont = [SYMBOL, 16],color=[red,green],tickmarks=[9,12],thickness=2):
P2:=implicitplot([sigma=curve2,sigma=curve1],sigma = -1.2 .. 2, A = A0 .. 2,  numpoints = 20000, axes = box, axesfont=[SYMBOL, 14],labels = [sigma, A], labelfont = [SYMBOL, 16],color=[red,blue],tickmarks=[9,12],thickness=2):
display(P1,P2);
                       

@Ronan 

fact := proc (expr::equation) options operator, arrow; frontend(gcd, [numer(lhs(expr)), numer(rhs(expr))])/frontend(gcd, [denom(lhs(expr)), denom(rhs(expr))]) end proc;
 
eq2 := (F+G)*R*sin(x)/sqrt(W-1) = (1/2)*M*R^2*tan(x)*(a/(R*sqrt(W-1)));
eq22:=convert(eq2, sincos);
 fact(eq22); eq22/fact(eq22);

@tdavid  If  rho(phi)  is known, then you must represent this function here. For this you can (after saving your document first) upload it here, using the bold green up-arrow in the MaplePrimes editor.

For example if a curve is an ellips  x^2/5^2+y^2/3^2=1  then there is no problems:

rho:=t->1.8/(1+0.8*cos(t));  # Polar equation of this ellips is  r=rho(t)  (the pole in the right focus)
L:=phi->Int(sqrt(rho(t)^2+diff(rho(t),t)^2), t=0..phi);
plot(L(phi), phi=0..2*Pi);
evalf(L(4));
evalf(L(5));
phi=fsolve(L(phi)=50);

@tzeng  Below is the corrected code:

restart;
Mphi := f->-I*(diff(f, phi));
Pr := f->-I*(diff(f, r));
rplus := f->r*exp(I*phi)*f;
rminus := f->r*exp(-I*phi)*f;
Pplus := f->exp(I*phi)*(Pr(f)+I*Mphi(f)/r);  
Pminus := f->exp(-I*phi)*(Pr(f)-I*Mphi(f)/r);
Rpp := f->(Pplus(f)+I*rplus(f));
Rmp := f->(Pplus(f)-I*rplus(f));
Rpm := f->(Pminus(f)+I*rminus(f));
mm := f->(Pminus(f)-I*rminus(f));
Psi00 := sqrt(1/Pi)*exp(-(1/2)*r^2);                        
Rpm(Rpp(Psi00));                
simplify(%);   
 

@Earl  Wonderful code! Thanks for this.