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

Usually, select or map produce shorter and faster solutions than a loop.

select(`>`, L, 0) /~ 2;

DD:=(n,alpha) -> Matrix(n+1, (i,k) -> sum(omega(i-1,j,k-1), j=ceil(alpha)..n)):
DD(2,alpha);

It remains to define your function omega:

omega := (i,j,l) -> ...

Yes, in Maple 2019 it is recommended to use

g := x-> local j; sum(..., j=1..5);