You want to simplify a polynomial f in n (=7) indeterminates.
Maple has simplify(f, size)  for this (size is by default).
Of course the size is not necessarily minimal, and AFAIK such algorithm does not exist.
But if you know/hope/suspect a special convenient form for f, tell us about it.

restart;
f:=
arccos((-sin(theta3)*cos(phi3)*cos(theta2)+cos(theta3)*sin(theta2))/((
(1-cos(theta3)^2)*cos(theta2)^2+cos(theta3)^2-1)*cos(phi3)^2-2*sin(theta3)*cos(
theta2)*cos(theta3)*sin(theta2)*cos(phi3)-cos(theta2)^2*cos(theta3)^2+1)^(1/2))
*sin(theta3):

evalf[15](Int(Re(f), [theta2=0..Pi,theta3=0..Pi,phi3=0..2*Pi], method=_CubaCuhre, epsilon=1e-6));
                        62.0124733278130
evalf[15](Int(Im(f), [theta2=0..Pi,theta3=0..Pi,phi3=0..2*Pi], method=_CubaCuhre, epsilon=1e-6));
                       Float(undefined)

See here.
It will be easy to use Explore and set your prefered sliders.

It is a bug in the symbolic part. The integral is simplified to BesselI for noninteger indices.

restart;
F := Int(exp(z*cos(theta))*cos(v*theta)/Pi, theta = 0 .. Pi):
ex := eval(F, [z=2, v=6/5]):
value(ex);     # wrong
    BesselI(6/5, 2)
evalf(%) = evalf(ex);
                   1.389981292 = 1.377481960
# FunctionAdvisor(BesselI, integral_form);

P.S. If you are intersted only in numeric results, just use Int instead of int and it will be OK.

It is used the fact that A,B are similar iff the characteristic matices (A-x*I, B-x*I) have the same Smith form.
So, the eigenvalues are not needed, and the computations are in Q if A,B are over Q.

A:=n ->
Matrix( n+1, (i,j)->
`if`(i<j,0,
  `if`(i=1,1,
  (i-1)*(-1)^(i-j)*2^(2*j-2)*(i+j-3)! / ( (2*j-2)!*(i-j)! )
))):

A(5);

You must use

PA[i] := eval(x[i + 1](t), sol(t[i]));

(note the i+1 because  x[0] is not defined) and similar for the rest of the loops.

Also, dt := 0.002;  because for 0.02 you will have singularities.

It is recommended to change tt[i] instead of t[i], because t is already in use.

For this series (with a numeric n) there are some subtle math aspects.

1. The series may not exist if a<n [e.g. if a is noninteger]
2. For a>=n the expansion is 0  i.e. 0 + O(...)
but unfortunately Maple is not able to find it because for example

diff(sin(x)^a, x);

and Maple tries x=0 instead of limit(%, x=0)

You may of course expand around another point, e.g.

series(sin(x)^a, x=Pi/2, 4);

The labels (A,B)  persist until the output is deleted. Even if the vector <1,2,3> is changed.

ws.mw

It is of course possible, but I don't think it's a good idea.

You may obtain animations using Explore.
Here is a simple example, similar to your problem.

restart;
A:=proc(c)
local p1,p2,d,x,t ;
p1:=plot(1 + c*sin(t), t=0..2*Pi, coords=polar);
d:=fsolve(x^5+x = c^5+1, x);
p2:=plot(2*x+d, x=-2..1);
plots:-display(p1,p2, view=[-2.5 .. 2.5, -1..4]) 
end proc:

Explore(A(c), c=0..2.6, animate);

BISSEC := proc (PP, UU, VV) 
uses LinearAlgebra,plots;
local P:=PP,U:=UU, V:=VV,
a, b, eq1, eq2, M1, M2, t, PU, PV, bissec1, bissec2; 
a := (P-U)/Norm(P-U, 2)+(P-V)/Norm(P-V, 2); M1 := P+a*t; b := (P-U)/Norm(P-U, 2)-(P-V)/Norm(P-V, 2); M2 := P+b*t; 
eq1 := op(eliminate({x = M1[1], y = M1[2]}, t)); 
eq2 := op(eliminate({x = M2[1], y = M2[2]}, t)); 
P := convert(P, list); U := convert(U, list); 
V := convert(V, list);
PU := plot([P, U]); PV := plot([P, V]); bissec1 := implicitplot(op(eq1[2]), x = 0 .. 5, y = 0 .. 10, color = red); bissec2 := implicitplot(op(eq2[2]), x = 0 .. 5, y = 0 .. 10, color = green); display([bissec1, bissec2, PU, PV], scaling = constrained) end proc:

A := `<,>`(4, 8): B := `<,>`(4, 2):; C := `<,>`(1, 4):
#debug(BISSEC);
BISSEC(A, B, C);

P:=n -> binomial(n-1,3) - binomial(floor(n/2),3);