Try this:

contourplot(psi(x, eta), x = 0 .. 5, eta = 0 .. 10, filledregions=true, coloring=["Blue", "Red"]); 
 

Another way is to use  Student:-Calculus1:-Roots  command. It finds symbolically all the real roots of some equation (your system is easily reduced to a single equation). Simultaneously some simplification of the found roots is made:

restart;    
f1:=x^2+y^2=1:    
f2:=y=x^3:    
Student:-Calculus1:-Roots(subs(f2,f1)):
seq([x,eval(y,f2)], x=%);  # All the real roots symbolically
evalf([%])[];  # All the real roots numerically
                  
 

Download Roots.mw

1. If you need to build the surface  f1  for a fixed point in time (for example for  x1=30 ), then you can do it this way:

restart;
x1:=30; 
f1:=11.921131-0.140022*x2+0.001034*x3+0.000836*x1*x1+0.001003*x2*x1-0.000052693*x2*x2+0.000177*x3*x1+0.000157*x3*x2-0.000011584*x3*x3;
plot3d(f1, x2=0..100, x3=400..800, axes=boxed);

2. If you want to investigate how this surface changes with changing time, then it is convenient to use  Explore  command:

restart;
f1:=11.921131-0.140022*x2+0.001034*x3+0.000836*x1*x1+0.001003*x2*x1-0.000052693*x2*x2+0.000177*x3*x1+0.000157*x3*x2-0.000011584*x3*x3:
Explore(plot3d(f1, x2=0..100, x3=400..800, axes=boxed), x1=15..60);

Moving the slider by the mouse, you can observe the change of the surface with the time change.

Probably means the plotting of three arrows in different colors. When you try to build this by one command, you actually get a similar error.

Here is an workaround:

with(plots):
Vectors:=[[1,0,0],[0,1,0],[0,0,1]]:
Colors:=[red,blue,green]:
display(seq(arrow(Vectors[i], shape = cylindrical_arrow, length = 1, color=Colors[i]), i=1..3));
                  

# An error below
arrow(Vectors, shape = cylindrical_arrow, length = 1, color=Colors); 
 

Use  Physics  package for this:

restart;
with(Physics):
Setup(noncommutativeprefix={A,B});
expand((A+B)^3);
 

It is not at all necessary to remember the Maple commands for calculating the norm and calling the corresponding packages. The main thing is to know how to introduce a vector and know the formula for the norm:

restart;
A:=<4, 7>;  # The shortest way to specify a column 2D-vector
sqrt(A[1]^2+A[2]^2);  # The calculation of the Euclidean norm symbolically
evalf(%);  # The same but numerically

If to solve as proposed in another answer, then we unjustifiably narrow the domain of the expression  A . In fact,  the identity  A = -b/a  is true for any real numbers (of course  a<>0). To work with a real  x  in the real domain, it's better to write surd(x, 3)  instead of  x^(1/3) . See the examples:
(-8)^(1/3);
simplify((-8)^(1/3));
surd(-8, 3);
                          

Using  surd  command, the initial expression  A  can be simplified in 2 steps:

restart;
A := surd(a^6,3)*surd(-b^3,3)/a^3;
simplify(A) assuming real, a>0;
simplify(A) assuming real, a<0;
                           

See help on  surd  command for details.

restart;
eqq:= k[t]*(`ℓ`^2)*(diff(q[3](tau), tau, tau)+(5*alpha-sigma+2*theta+1)*q[3](tau)+(-4*alpha+sigma-theta)*q[2](tau)+q[1](tau)*alpha) = -(sqrt(m*(1/k[t]))*`ℓ`*k[t]*`Δθ`*(q[3](tau)-q[2](tau))*sin(sqrt(Lambda*k[t]*(1/m))*sqrt(m*(1/k[t]))*tau)+2*xi*sqrt(lambda*k[t]*m)*(diff(q[3](tau), tau)))*`ℓ`*(1/sqrt(m*(1/k[t]))):
vars:= [diff(q[1](tau),tau$2),diff(q[2](tau),tau$2),diff(q[3](tau),tau$2),diff(q[1](tau),tau),diff(q[2](tau),tau),diff(q[3](tau),tau),q[1](tau),q[2](tau),q[3](tau)];
expand(simplify((eqq-rhs(eqq))/(k[t]*`ℓ`^2))) assuming positive;
collect(%, vars);

Edit.

In Maple 2018 use empty symbol  ``  for this. In this case, Maple regards  ``(1+A*cos(psi))
as a single expression:

restart;
signum(``(1+A*cos(psi))-sqrt(1+sin(psi))) assuming ``(1+A*cos(psi))<sqrt(1+sin(psi));
                                                               -1

Replace  the line 

pp1:=display([pu1,pu2,pu3,pu4]); 

with

pp1:=plots:-display([pu1,pu2,pu3,pu4]);

It is shorter and faster to use  plot  command for this (without calling  plots  package):

L1:=[seq(j, j=1..10)];
L2:=[seq(j^2, j=1..10)];
plot(L1, L2, style=point, symbol=solidcircle, symbolsize=20);
# Or
plot(L2, L1, style=point, symbol=solidcircle, symbolsize=20);

 #  If you omit  style=point  option, then (by default) Maple draws a polyline that connects these points:

plot(L1, L2);
 

Your system has no solutions, since we have

dsolve({diff(x(t), t) = y(t), diff(y(t), t) = -x(t), x(0) = 0, y(0) = 1});
evalf([(D(sin))(8), (D(-cos))(8)]);
                                {x(t) = sin(t), y(t) = cos(t)}
                          [-0.1455000338, -0.9893582466]
 

You must impose the condition of touching the circles  C4  and  C2 : the distance between the centers is equal to the sum of the radii. Here is the solution:

restart; with(geometry):

x0 := (3+1/2)/2;
x1 := x0; y1 := x0 + (1+1/4);
point(P1, x1, y1); ## center of the circle
circle(C2, [P1, x0]);

x3 := (3+1/2)/2-(3/4+20/1000)/2-1/2; y3 := 1+1/4; evalf([x3, y3]);

point(P2, x3, 0); point(P3, x3, y3);
line(L1, [P2, P3]);
intersection(I1, L1, C2);
for s in I1 do print(evalf(coordinates(s))) end do;
tr := 1/8;  ## tool radius, radius of C4

x4 := x3 - tr;

point(P4, x4, y4);
circle(C4, [P4, tr]);
Equation(C4,[x,y]);
distance(P4,P1);
solve({%=tr+x1,y4<y1});
point(P4,x4,eval(y4,%));
circle(C4, [P4, tr]);
intersection(I2, C2, C4);  # The point of touching the circles  C4  and  C2
detail(I2);
draw([L1,C2,C4], view=[0..4,0..5], axes=normal);

I usually prefer to use lists or vectors; over them it is more convenient to make various arithmetic manipulations:

N:=10:
Y := [seq(i^2, i=1 .. N)];   # This is the list
DY := Y[3 .. N] - Y[1 .. N-2];  # This list also

# or the same, but slightly longer:

DY := [seq(Y[n+1]-Y[n-1], n=2..N-1)];

sqrt((360-alpha[1]-3*alpha[2])/``(720/Pi^2*(cos(Pi*alpha[1]/180)+cos(Pi*alpha[2]/180))^2)-1);

#  or

sqrt((360-alpha[1]-3*alpha[2])/(``(720/Pi^2)*(cos(Pi*alpha[1]/180)+cos(Pi*alpha[2]/180))^2)-1);