with(LinearAlgebra):with(plots):with(plottools):
A1:=<0,0,1>: A2:=<0,3/5,4/5>: A3:=<4/5,3/5,0>: A4:=<5/13,0,12/13>:
seg:= (a,b)->spacecurve(Normalize(a+t*(b-a),2),t=0..1, thickness=3, color=red):
display(sphere([0,0,0],1),seg(A1,A2),seg(A2,A3),seg(A3,A4),seg(A1,A4));

You can use e.g.

binaryvariables = indets(L, specindex(integer,x));

where L is a list containing the constraints (and the function).
Or, you may generate their sequence using seq.
 

f:=(x,y) ->exp(-x^2-y^2):
L:=4: n:=64:
r1:=rand(0..1.0): r:=rand(-1..1.0):
a:=Vector(n,L*r): b:=Vector(n,L*r):
c:=Vector(n,0.3*r):
da:=Vector(n, i->0.5*''r1()''):
P:=proc(aa)
  global n,a,b,c,da,h,L;
  a:=a+eval~(da);
  map[inplace](frem,a,L);
  h:=add(c[i]*f(x+a[i], y+b[i]),i=1..n):
  plot3d(h, x=-L..L,y=-L..L, color="Aquamarine", scaling=constrained, view=-2..2, axes=none, style=surface, lightmodel=light1);
end:
plots:-display(seq(P(aa),aa=1..50),insequence=true);
## or ...
Explore(P(aa), aa=0..1, animate,loop,size=[800,800]);

Download req-int.mw

restart;
f:=r^2 *cos(theta)+r*sin(theta):
r1:= 0.3+0.1*cos(theta):
r2:= 0.5+0.1*cos(theta):
plots:-densityplot(f, theta=0..2*Pi, r=r1 .. r2,  
coords=polar, colorstyle = HUE, style = patchnogrid);

If you need labels:

plots:-display(%, labels=["x","y"]);

c:=Matrix([    # result = 22.613, [4, 2, 5, 3, 6, 1, 4]
[0., 3.60555127546399, 2.00000000000000, 8.06225774829855, 5.09901951359278, 1.41421356237310],
[3.60555127546399, 0., 2.23606797749979, 6.32455532033676, 2.23606797749979, 3.60555127546399],
[2.00000000000000, 2.23606797749979, 0., 8.06225774829855, 3.16227766016838, 1.41421356237310],
[8.06225774829855, 6.32455532033676, 8.06225774829855, 0., 8.06225774829855, 9.00000000000000],
[5.09901951359278, 2.23606797749979, 3.16227766016838, 8.06225774829855, 0., 4.47213595499958],
[1.41421356237310, 3.60555127546399, 1.41421356237310, 9.00000000000000, 4.47213595499958, 0.]]):

n := sqrt(numelems(c)):  N:={seq(1..n)}:
z:=add(add( c[i,j]*x[i,j],j=N minus {i}),i=N): 
conx:=seq( add(x[i,j], i=N minus {j})=1,j=N),  seq( add(x[i,j], j=N minus {i})=1,i=N):
conu:= seq(seq(u[i]-u[j]+n*x[i,j] <= n-1, i=N minus {1,j}), j=N minus {1}):

r := Optimization[LPSolve](z, {conx, conu}, #integervariables={seq(u[i],i=N minus {1})},
                           binaryvariables={ seq(seq(x[i,j],i=N minus {j}),j=N)} );

X:=eval(Matrix(n, symbol=x),{r[2][], seq(x[i,i]=0,i=1..n)}):
ord:=1: k:=1: to n do k:=max[index](X[k]); ord:=ord,k od:'order'=ord;

r := [22.6135858310497, [u[2] = -.999999998967403, u[3] = -2.99999997900383, u[4] = -8.88178419700125*10^(-16), u[5] = -2.00000000722817, u[6] = -3.99999998106902, x[1, 2] = 0, x[1, 3] = 0, x[1, 4] = 0, x[1, 5] = 0, x[1, 6] = 1, x[2, 1] = 0, x[2, 3] = 0, x[2, 4] = 1, x[2, 5] = 0, x[2, 6] = 0, x[3, 1] = 0, x[3, 2] = 0, x[3, 4] = 0, x[3, 5] = 1, x[3, 6] = 0, x[4, 1] = 1, x[4, 2] = 0, x[4, 3] = 0, x[4, 5] = 0, x[4, 6] = 0, x[5, 1] = 0, x[5, 2] = 1, x[5, 3] = 0, x[5, 4] = 0, x[5, 6] = 0, x[6, 1] = 0, x[6, 2] = 0, x[6, 3] = 1, x[6, 4] = 0, x[6, 5] = 0]]

order = (1, 6, 3, 5, 2, 4, 1)

If you need to read the same data several times you could do this:

1. Import once using ImportMatrix into A
2.    save A, "hugeA.m";
3.  Whenever you need A, use
     read "hugeA.m";
     It will be executed much faster (being in internal format .m).

The result is wrong, the implicitplot algorithm does not work well here.




Probably some "clever" options could be added, but always a better idea is to try a parametric representation of the curve:

restart;
f:=(1/2)*x*tan(log((1/5)*(2*(x^2+y^2)))) = y:
g:=simplify(eval(f, [x=r*cos(t), y=r*sin(t)])) assuming r>0:
u:=solve(g,t);

plot( [[r*cos(u),r*sin(u), r=0..5],[-r*cos(u),-r*sin(u), r=0..5]]);

The idea is similar to Kitonum's but has some advantages.
Note that some details will be omitted.

As it was observed,

f:=a^2+b^2+c^2+a*b+a*c+b*c-(a+b-c)*sqrt(2*a*b+a*c+b*c)-(a+c-b)*sqrt(a*b+2*a*c+b*c)-(b+c-a)*sqrt(a*b+a*c+2*b*c);
is positively homogeneous and symmetric, so we may suppose:
-1 <= a, b, c <= 1.

It is easy to see that c=0 implies a=b=0.
It will be enough to consoder c=1 and c=-1.

1.     c=1

f1:= eval(f, c=1);
plots:-implicitplot(f1, a=-1..1,b=-1..1, grid=[1000,1000]);

So, we have two branches. But they are symmetric because f(a,b,1)=f(b,a,1).
So we may consider only a>0, b<0.
Actually a will be in a1..a2  where a1, a2 could be computed [here at least of one sqrt must be 0].
So, the corresponding b will be h(a) = RootOf(f1, b=-1/5 .. -1/3).
[Note that h(a) can be expressed as a RootOf of a polynomial]

2.    c=-1
There are no solutions.

Conclusion:
The roots are:  a=k*u, b=k*h(u), c=k,  for all u in [a1,a2] and k>=0
and all permutations of these.

Edit.
 

a1:=-1+sqrt(1+RootOf(_Z^8-8*_Z^7-22*_Z^6+228*_Z^5+17*_Z^4-624*_Z^3-352*_Z^2+256*_Z+256, index = 2));

a2:=1;

h:= a -> 
RootOf(
a^16-8*a^15*b-4*a^14*b^2-48*a^13*b^3+154*a^12*b^4+40*a^11*b^5-304*a^10*b^6-56*a^9*b^7+451*a^8*b^8-56*a^7*b^9-304*a^6*b^10+40*a^5*b^11+154*a^4*b^12-48*a^3*b^13-4*a^2*b^14-8*a*b^15+b^16-8*a^15-8*a^14*b-192*a^13*b^2-56*a^12*b^3+96*a^11*b^4+456*a^10*b^5+296*a^9*b^6-560*a^8*b^7-560*a^7*b^8+296*a^6*b^9+456*a^5*b^10+96*a^4*b^11-56*a^3*b^12-192*a^2*b^13-8*a*b^14-8*b^15-4*a^14-192*a^13*b-480*a^12*b^2-1104*a^11*b^3+1564*a^10*b^4+3040*a^9*b^5-404*a^8*b^6-4048*a^7*b^7-404*a^6*b^8+3040*a^5*b^9+1564*a^4*b^10-1104*a^3*b^11-480*a^2*b^12-192*a*b^13-4*b^14-48*a^13-56*a^12*b-1104*a^11*b^2+720*a^10*b^3+4616*a^9*b^4+3464*a^8*b^5-4864*a^7*b^6-4864*a^6*b^7+3464*a^5*b^8+4616*a^4*b^9+720*a^3*b^10-1104*a^2*b^11-56*a*b^12-48*b^13+154*a^12+96*a^11*b+1564*a^10*b^2+4616*a^9*b^3+6008*a^8*b^4-1792*a^7*b^5-7416*a^6*b^6-1792*a^5*b^7+6008*a^4*b^8+4616*a^3*b^9+1564*a^2*b^10+96*a*b^11+154*b^12+40*a^11+456*a^10*b+3040*a^9*b^2+3464*a^8*b^3-1792*a^7*b^4-3896*a^6*b^5-3896*a^5*b^6-1792*a^4*b^7+3464*a^3*b^8+3040*a^2*b^9+456*a*b^10+40*b^11-304*a^10+296*a^9*b-404*a^8*b^2-4864*a^7*b^3-7416*a^6*b^4-3896*a^5*b^5-7416*a^4*b^6-4864*a^3*b^7-404*a^2*b^8+296*a*b^9-304*b^10-56*a^9-560*a^8*b-4048*a^7*b^2-4864*a^6*b^3-1792*a^5*b^4-1792*a^4*b^5-4864*a^3*b^6-4048*a^2*b^7-560*a*b^8-56*b^9+451*a^8-560*a^7*b-404*a^6*b^2+3464*a^5*b^3+6008*a^4*b^4+3464*a^3*b^5-404*a^2*b^6-560*a*b^7+451*b^8-56*a^7+296*a^6*b+3040*a^5*b^2+4616*a^4*b^3+4616*a^3*b^4+3040*a^2*b^5+296*a*b^6-56*b^7-304*a^6+456*a^5*b+1564*a^4*b^2+720*a^3*b^3+1564*a^2*b^4+456*a*b^5-304*b^6+40*a^5+96*a^4*b-1104*a^3*b^2-1104*a^2*b^3+96*a*b^4+40*b^5+154*a^4-56*a^3*b-480*a^2*b^2-56*a*b^3+154*b^4-48*a^3-192*a^2*b-192*a*b^2-48*b^3-4*a^2-8*a*b-4*b^2-8*a-8*b+1,
b, -3/10 .. -27/100)   :

To plot the solution use:

plot3d([ [k*u, k*h(u), k], [k*u, k, k*h(u)], [k*h(u), k*u,  k],
[k*h(u), k, k*u], [k, k*h(u), k*u], [k, k*u, k*h(u)] ],  u=a1..a2, k=0..10);

A strange fact (which I don't understand) is that in Maple there are both left and right versions of spherical coordinates!

The left one is for plot3d:
[u,v,w] -> [u*cos(v)*sin(w), u*sin(v)*sin(w), u*cos(w)]

and the right one is in VectorCalculus:
[u,v,w] -> [u*cos(w)*sin(v), u*sin(w)*sin(v), u*cos(v)].

Of course, as Carl has shown they can be redefined by the user, but it is confusing.

restart;
S:= rhs~(solve({-x-y+3*z =a, -x+2*y = b, 3*x-2*y-z = c}, [x,y,z])[]):
A:=eval(S, [a=1,b=0,c=0]):
B:=eval(S, [a=0,b=1,c=0]):
C:=eval(S, [a=0,b=0,c=1]):
OO:=[0,0,0]:
plots:-polygonplot3d([[OO,A,B], [OO,A,C], [OO,B,C]]);

In your example, you can do this only if the rank of the system is r=19 and the rank corresponding to the unknowns Z1,...,Z19 is also r=19.
In this case you can solve wrt Z1,...,Z19, e.g.
solve( sys, {seq(Z||i,i=1..r)} );