The  coeff  command does not work for extracting the coefficients of polynomials from several variables. A special procedure is required for this. The  coefff  procedure extracts the coefficient in front of a monomial  t of the polynomial  P  from the variables  T .
 

Download coefff.mw

Edit.

In the first equation there should be the multiplication sign after  x , otherwise everything that stands in parentheses (after  x)  disappears.
See:

x:=2:
x(h(t));
                                 2

TM:=proc(N)
if N<2 then error "Should be N>=2" fi;
Matrix(N+1, [[seq(p[i],i=1..N-1),o[3]], [o[1],seq(q[i],i=1..N-1),o[4]], [o[2],seq(r[i],i=1..N-1)]], shape = band[1,1], scan=band[1,1]);
end proc:

Example of use:

TM(5);
                                

To multiply the matrices, I replaced  &*  with the dot  . . Now everything works.

Help_(3)_new.mw

Edit.

for this. 

Example:

plot([cos(t), sin(t), t=-1..1], x=0..2, y=-2..2, color=red, scaling=constrained);   # Or

plot([cos(t), sin(t), t=-1..1], color=red, scaling=constrained, view=[0..2,-2..2]); 
 

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.