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);

Maple can not solve your equation  deq=3  for  x , since the explicit dependence of  y(x)  is unknown. Here is a workaround:

restart:
eq:=4*x^2+2*y(x)^2=32.5625:
deq:=solve(  diff(eq,x), diff(y(x),x) );
deq=3;
y(x)=solve(deq=3, y(x));
x=solve(subs(y(x)=y, deq=3), x);

RootOf(2*_Z+3*y(_Z))  is a formal representation of the root of the equation  2*_Z+3*y(_Z)=0  with respect to  _Z

I corrected all the syntax errors. We see that the integrand is actually a scalar function.  But Maple does not cope with this integral even for finite intervals. I think that the reason is in a large number of parameters (8 parameters).

You can try to calculate your integral by specifying the values of these parameters.

Edit.

Download Integration-Vec-Example_new.mw

restart; 
with(Physics):
Setup(mathematicalnotation = true);                
Setup(signature = `-+++`);                    
Coordinates(M = [t, rho, z, phi]);
ds2 := -exp(2*psi(rho, z))*dt^2+exp(2*gamma(rho, z)-2*psi(rho, z))*(drho^2+dz^2)+exp(-2*psi(rho, z))*(rho^2)*(dphi^2);  
Setup(metric = ds2);

Download subs.mw