plot( x[3]^5, caption = typeset("\n A plot of %1.", x[3]^5), captionfont=[times, 20] );

We use the classical method of finding critical points and their subsequent investigation with the help of the second derivative. Maple returns the results in thу range  -Pi .. Pi . If you want the results to be in the range  0..2*Pi , just add  2*Pi  to the corresponding variable:  

restart;
f:=(x,y)->1+8*cos(1/2*x-1/2*y)*cos(1/2*x)*cos(1/2*y);
[solve({D[1](f)(x,y)=0, D[2](f)(x,y)=0},explicit)];
S:=[seq([Student:-MultivariateCalculus:-SecondDerivativeTest(f(x,y),[x,y]=eval([x,y],%[i]))], i=1..nops(%))];
select(s->s[1]<>(LocalMin = []), S);

The final results:

with(plots):
K:=9;
deG:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*sin(theta(t))= 0;
deL:=diff(theta(t),t,t) + mu*diff(theta(t),t)+K*theta(t)= 0;
Iv:=theta(0)=0.75, D(theta)(0)=2.0;
dom1:=t=0..10;
soln1a:=dsolve({eval(deL,mu=0),Iv});
gr1a:=plot(eval(theta(t), soln1a), dom1, color=blue);
 

Your system, in addition to the unknown, contains several parameters. You can obtain explicit solutions depending on the values of these parameters as follows:

Download System_with_parameters.mw

x^2:
f:=unapply(%, x);
f(t);
f(3);

We can find all 3 solutions as follows:

restart;
z:=x+I*y:
evalc(abs(z)*(z-4-I)+2*I = (5-I)*z);
solve({Re(lhs(%))=Re(rhs(%)), Im(lhs(%))=Im(rhs(%))}) assuming real;
evalf(%);
                     

Edit. I do not know why  solve  (as OP did) fails with this example. Probably this should be seen as a bug.

Maple plots everything correctly. What you call a solid disk is just a plot of the zero function, because the square root of a negative number is a purely imaginary number and its real part is 0 , for example

sqrt(-2);
Re(%);
                                   I*sqrt(2)
                                       0            

In addition, replace  c=-1..1  in your code with  c=-1...1. , otherwise the Maple will plot graphics for only three integer values  c=-1, 0, 1                

I do not know how the Groebner bases are related to your problem, but if you just need to find the minimum number of buses that satisfy the conditions, then you can do so:

Optimization:-LPSolve(add(x || i, i = 1 .. 30), B30[1 .. 30], assume = binary);

The result:

 [10, [x1 = 0, x10 = 1, x11 = 1, x12 = 1, x13 = 0, x14 = 0, x15 = 0, x16 = 0, x17 = 0, x18 = 0, x19 = 1, x2 = 1, x20 = 0, x21 = 0, x22 = 0, x23 = 0, x24 = 1, x25 = 0, x26 = 1, x27 = 0, x28 = 0, x29 = 1, x3 = 1, x30 = 0, x4 = 0, x5 = 0, x6 = 1, x7 = 0, x8 = 0, x9 = 0]]
 

Here are 2 procedures that do the same things as vv's ones, but with open source:

Two_procedures.mw

int(rationalize(1/(1+sqrt(x))), x=0..1);
                                                2-2*ln(2)

Explanation: rationalize command allows you to get rid of the radicals in the denominator, for example:

rationalize(1/(sqrt(2)-1));
                                                1+sqrt(2)

When simplifying expressions with radicals, it is often useful to combine  rationalize  with  simplify  or  expand :

simplify(rationalize(2/(sqrt(2)-1)+2/(sqrt(2)-2)));
expand(rationalize(2/(sqrt(2)-1)+2/(sqrt(2)-2)));
                                                  sqrt(2)
                                                  sqrt(2)

​​​​​​​Edit.

f:= x-> (sign(4*x^2-4)*abs(4*x^2-4)^(1/5))^4;
evalf(f(0));

Use  plots:-display  command for this:

plots:-display(rect[i], i=1..21);

Edit: of course above should be
plots:-display(seq(rect[i], i=1..21));

Thanks to Carl for the amendment.
 

You can treat this as a curve in 3d.

Example:

restart;
eq1:=sin(z):
eq2:=cos(z):
plots:-spacecurve([eq1,eq2,z], z=0..4*Pi, color=red, thickness=3);

Everything works with assuming:

simplify(A) assuming real;
simplify(B) assuming real;
simplify(A) assuming complex;
simplify(B) assuming complex;
combine(A);
combine(B);
combine(A*B);

                                                          1
                                                          1
                                                          1
                                                          1
                                                          1
                                                          1
                                                          1

Or

simplify(numer(A))/denom(A);

and so on ...

Several other ways:

is(A=1);
is(B=1);
is(A*B=1);
combine(convert(A, tan));
combine(convert(B, tan));
combine(convert(A*B, tan));

Edit.

I swapped  seq  and  plottools:-arrow  and added a few required parameters for  plottools:-arrow  command.

plot_arrows_sequence_new.mw

Addition.  [seq] means that the result of applying  seq  command will be surrounded by square brackets, i.e. we get a list.

Example:

seq(k^2, k=1..10);
[seq](k^2, k=1..10);