But your code works.
The only thing I did was to copy it in a new worksheet (1d math, worksheet mode) and fill some ";",
because I never debug in document mode.

restart;
with(plots):
with(plottools):
Cylinder := proc (a::list, b::list, dia := .1, colr := "Blue"): 
plots:-arrow(a, b, width = dia, shape = cylindrical_arrow, head_length = 0., head_width = 0.1e-1, colour = colr) 
end proc:
Body2 := display([point([1, -.2, 0], colour = red, symbolsize = 50, symbol = solidsphere), point([-1, -.2, 0], colour = red, symbolsize = 50, symbol = solidsphere), point([0, 1, 0], colour = purple, symbolsize = 50, symbol = solidsphere), Cylinder([0, 0, 0], [1, -.2, 0], 0.75e-1, "Red"), Cylinder([0, 0, 0], [1, -.2, -0.1e-1], 0.75e-1, "Nautical GrayViolet"), Cylinder([0, 0, 0], [-1, -.2, 0], 0.75e-1, "Red"), Cylinder([0, 0, 0], [-1, -.2, -0.1e-1], 0.75e-1, "Nautical GrayViolet"), Cylinder([0, 0, 0], [0, 1, 0], 0.65e-1, "Green"), Cylinder([0, 0, 0], [0, 1, -0.1e-1], 0.65e-1, "Nautical GrayViolet"), line([0, 0, 0], [.5, 0, 0], colour = red, thickness = 4), line([0, 0, 0], [0, .5, 0], colour = yellow, thickness = 4), line([0, 0, 0], [0, 0, .5], colour = blue, thickness = 4), plots:-textplot3d([.5, 0.5e-1, 0, "X"]), plots:-textplot3d([.1, .5, 0, "Y"]), plots:-textplot3d([0.5e-1, 0, .5, "Z"])], scaling = constrained);
w1 := .35; w2 := .45; w3 := .15;
animate( rotate, [Body2,w1*t,w2*t,w3*t], t=0..8*Pi,frames=100 );

For rotate alone, you must use e.g.

t:=Pi;
display(rotate(Body2, -w1*t, -w2*t, -w3*t));

rota.mw

Are you sure that the ODE has a solution? Explore suggests that it has not.

restart;
Digits := 15:
alpha := evalf(2*Pi*(1/180)); EP := .2; lambda := .1; HA := 5; RE := 20;
ODE := (1+EP)*(((D@@3)(F))(r)+4*alpha^2*(D(F))(r))+2*alpha*RE*F(r)*(D(F))(r)-HA*alpha^2*(D(F))(r)-3*EP*lambda*((1/2)*(D(F))(r)^2*((D@@3)(F))(r)+(D(F))(r)*((D@@2)(F))(r)^2)/alpha^2-EP*lambda*(72*F(r)^2*(D(F))(r)+2*(D(F))(r)^3+32*F(r)*(D(F))(r)*((D@@2)(F))(r)+2*F(r)^2*((D@@3)(F))(r)) = 0:
#BCSforNum := F(0) = 1, (D(F))(0) = 0, F(1) = 0;
BCS :=        F(0) = 1, (D(F))(0) = 0, (D@@2)(F)(0) = a:
sol := dsolve({ODE, BCS}, numeric, parameters=[a], 'output'='listprocedure'):
P:=proc(A)
  sol('parameters'=['a'=A]); sol;
  plots:-odeplot(sol, [r, F(r)], r = 0 .. 1) 
end proc:
Explore(P(A),'parameters'=[A=-100. ..  100.],initialvalues=[A=0] );

You don't need Maple: the function is odd, so, the integral is 0.

The equation is not "simple".
1. It is not given explicitely y'(x) = f(x, y(x))
2. It has 2 conditions instead of one.

dsolve finds symbolically the correct solution:

dsolve(dsys);
     y(x)=1
because by chance it has one (even if two condition were used).
For a numerical solution you would need to give the explicit form and only one condition.
[Note that even in this case the equation  y' = sqrt(y-1) is not simple  because RHS is not Lipschitz near y0=1].

This is normal, sqrt is not differentiable at 0.

If you want an expansion (non-Taylor), you may use:

eval(mtaylor(tanh(B*J*Q*z+B*Jo*m), [Q, m], 9) , Q=sqrt(q));

See  ?sum,details  for the definition of sum when lowerlimit > upperlimit.

In your case
sum(-f(k+1,x), k = 0 .. -2) ;
    f(0,x)

So, your sum equals GAMMA(0,x)  =  Ei(1,x).

If you have a finite sum such as  s(x) = a[0] + a[1]*x + ... + a[9]*x^9
it is not possible to find a (unique) function whose Taylor (or other) series starts  with s(x).
[A similar problem is to find "the" irrational number for which the first decimals are 1.234]

Using your random generators,  the equality fails in about 2% of the tests.

restart;
with(RandomTools):
NN := 1000:
err := 0: 
to NN do 
d := Generate(integer(range = 1 .. 120)); 
T := Generate(integer(range = 1 .. d)); 
B := Generate(integer(range = 1 .. d)); 
Y := Generate(integer(range = 0 .. B)); 
Pr := Generate(integer(range = 0 .. d-T)); 
py := Generate(integer(range = 0 .. min(Y, d-T, Pr))); 
e1 := binomial(Y, py)*(binomial(d-Y, Pr-py)*binomial(d-Pr, T)-binomial(d-B, T)*binomial(d-T-Y, Pr-py))/(binomial(d-T, Pr)*(binomial(d, T)-binomial(d-B, T))); 
e2 := add((binomial(Y, n)*binomial(d-Y, T-n)/(binomial(d, T)-binomial(d-B, T))+n-max(n, binomial(d-B, T)/(binomial(d, T)-binomial(d-B, T))))*binomial(Y-n, py)*binomial(d-T-Y+n, Pr-py)/binomial(d-T, Pr), n = 0 .. Y); 

if e1 <> e2 then err := err+1; print(e1 <> e2, 'd'=d, 'T'=T, 'B'=B, 'Y'=Y, 'Pr'=Pr, 'pr'=py) end if end do: 
fails=err,'NN'=NN;

It is not a good idea to use floats in symbolic computations (if you can avoid them).

Replacing 0.2 by 1/5 (or convert( ...,rational) ), everything is OK.

Alternatively, increase Digits.

with(plots):
with(plottools):
Body := display([
point([1, -.2, 0], colour = red, symbolsize = 50, symbol = solidsphere), 
point([-1, -.2, 0], colour = red, symbolsize = 50, symbol = solidsphere), 
point([0, 1, 0], colour = purple, symbolsize = 50, symbol = solidsphere), 
line([0, 0, 0], [1, -.2, 0]), 
line([0, 0, 0], [-1, -.2, 0]), 
line([0, 0, 0], [0, 1, 0]), line([0, 0, 0], [.5, 0, 0], colour = red, thickness = 4), 
line([0, 0, 0], [0, .5, 0], colour = yellow, thickness = 4), 
line([0, 0, 0], [0, 0, .5], colour = blue, thickness = 4)], scaling = constrained):
w1 := .35: w2 := .45: w3 := .15:
animate( rotate, [Body,w1*t,w2*t,w3*t], t=0..8*Pi,frames=100 );

The cure is simple: define ff1 etc via unapply, e,g.

ff1 := unapply( display([line([0, 0, 0], [d[1][1], d[1][2], d[1][3]], colour = red, thickness = 4)]),  t);

Your code works. The only problem is that a,b,c,x3  are complex (in general) so in printf you must replace for example  %2.8f  with %2.8Zf  (and without spaces!).

You must assign ranges for all the variables (or for none).
You may set of course e.g.  A0 = -infinity .. infinity  if you have not a better range.

It is unlikely that the print occurs only at the end, unless a procedure consumes almost all of the running time.
Maybe you must scroll the page while running in order to see the messages.

Anyway, you can set e.g.
printlevel:=10;

before running the procedure and see exactly what it is doing. If necessary increase slightly the printlevel value.

SmithForm does not accept polynomials with complex coefficients.
But there exists LinearAlgebra:-Generic:-SmithForm which works in any Euclidean domain. In order to use it in C[z], some simple procedures must be provided by the user (see the help page). For example:

C[EuclideanNorm] := a -> degree(a,z):
C[Gcdex] := (a,b,s,t)->gcdex(a,b,z,s,t):
...

After that, you can use:
LinearAlgebra:-Generic:-SmithForm[C](A, output=['S','U','V']):