@Earl   See  emptysymbol  help page for this.

@vv  In fact, you have built a very thin ellipse. Once we realize the meaning of the problem, we can make it easier:
 plot([[1,1], [5,4]]);

@Rouben Rostamian  Thanks for pointing this out. Unfortunately, I do not know  output=...  option well and have never used it. However, I have never met an example where this option would be necessary. I think it is good that many tasks in Maple can be solved in different ways. We see that, for example, plots:-spacecurve command may well be replaced by  plots:-pointplot3d (with  style=line  option) or  plottools:-curve   commands.

@Carl Love  Thanks for the detailed explanation.

@jfreeman1212 

It's easy - you just need to solve the system:

solve({(x^2+y^2+2*y)^2=4*(x^2+y^2), y=x}, {x,y}, explicit);
           {x = 0, y = 0}, {x = sqrt(2)-1, y = sqrt(2)-1}, {x = -1-sqrt(2), y = -1-sqrt(2)}
                     

@Carl Love  Your procedure does not do what OP wants (look at  ListY  in the original post). In addition, it does not work in earlier versions than Maple 2018 (because of until).

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