"for loop" is not the best way to solve the problem. Here are 3 other ways:

restart;
A:=[1,2,3,4,5,6,7,8,9,10]:
B:=[10,20,30,40,50,60,70,80,90,100]:
n:=nops(A):
E:=5:
C:=a*E/4*Pi*b:

# Solution 1
DD:=[seq(eval(C,[a=A[i],b=B[i]]), i=1..n)];

# Solution 2
DD:=zip((u,v)->eval(C,[a=u,b=v]), A,B);

# Solution 3 
C:=(a,b)->a*E/4*Pi*b:
DD:=zip(C,A,B);

The number 3.14 ... in Maple is encoded as  Pi  not  pi. I also replaced D with DD since D cannot be a name (D is the differentiation operator in Maple).
We see that the shortest way is to use the  zip  command and write  C  as a function of 2 variables.

Download MaplePrimes_Example_anim.mw

restart;
Expr:=n->`*`(seq(seq(alpha[j]-alpha[i], j=i+1..n), i=0..n-1));

Example of use:
Expr(3);
                

 tomleslie's method does not work correctly for  t<0 . Below is a fixed and improved version. We use discont option to remove vertical lines where the function is discontinuous. These places can be identified through a dotted line (second plot). The values of the function at points multiple of  2*Pi  can be shown with small solid circles (third plot) :

Download periodic_function.mw

First, we find at what height ( z-value) these surfaces intersect, and then plot them using parametric equations and polar coordinates on the plane:

# z=sqrt(x^2+y^2)
solve({12-z^2=z,z>0});
A:=plot3d([r*cos(t),r*sin(t),12-r^2],r=0..3,t=0..2*Pi,color=khaki):
B:=plot3d([r*cos(t),r*sin(t),r],r=0..8,t=0..2*Pi,transparency=0.5):
plots:-display(A,B, scaling=constrained, axes=normal, labels=[x,y,``],labelfont=[times,14]);

b := 3*(4*x+1)/(diff(f(y),y)*(3*x+1)):
A:=(4*x + 1)/diff(f(y),y) = a:
subs(diff(f(y),y)=solve(A,diff(f(y),y)), b);

Execute

showstat(LinearAlgebra:-IsSimilar);

Example:

restart;
N:=3:
z:=x+I*y:
evalc(z^N);
                                 x^3-3*x*y^2+I*(3*x^2*y-y^3)

Examples.

Code in 1d math:
A:=<1,2; -1,5>;
LinearAlgebra:-Determinant(A);

Code in 2d math:
A:=<1,2; -1,5>;
|A|;

You probably have 2d math by default. Therefore, use the second example.

Edit. If you are not allowed to use a Maple procedure, here is a user-defined recursive procedure that evaluates the determinant of a matrix:

MyDet:=proc(A::Matrix)
local m,n;
m,n:=op(1,A);
if m<>n then error "The Matrix should be square" fi;
if m=1 then return A[1,1] else
add(A[1,k]*(-1)^(1+k)*MyDet(A[2..m,[seq(`if`(j<>k,j,NULL), j=1..n)]]),k=1..n) fi;
end proc:

Example of use:
A:=<1,2,3;4,5,6;7,8,10>;
MyDet(A);
 

You can calculate this integral numerically using a simple procedure as follows:

Int1 := (b,k,r,R)->evalf(Int(exp(-z*(R^2*k^2 - b^2*z)/(R*b))/(z*HeunB(0, k^2*R^2/(b*sqrt(R*b)), R^3*k^4/(4*b^3), 0, -sqrt(R*b)*z/R)^2), z = R .. r));

Example of use:

Int1(0.5, 1, 2, 0.9);
                                                1.462594758

Addition. Below I made an animation of this travel (continuation of the above code). The first frame shows the initially given points  A  and  C  and a green vector that determines the direction of movement:

Back:=display(textplot3d([cA[],"A"], font=[times,18], align=right),textplot3d([cC[],"C"], font=[times,18], align=right),arrow([1,2,4],v,color=green),pointplot3d([cA,cC],symbol=solidsphere,color=red,symbolsize=17)):
Anim1:=animate((display@plottools:-line),[cA,cA+~t*~(cB-cA), color=red,thickness=3],t=0..1,background=Back,paraminfo=false):
Anim2:=animate(display,[textplot3d([cB[],"B"], font=[times,18], align=right)],a=0..1,frames=5,background=display(Back,pointplot3d(cB,symbol=solidsphere,color=red,symbolsize=17)),paraminfo=false):
Anim3:=animate((display@plottools:-line),[cB,cB+~t*~(cC-cB), color=red,thickness=3],t=0..1,background=display(Back,pointplot3d(cB,symbol=solidsphere,color=red,symbolsize=17)),paraminfo=false):
display([Anim1, display(op([1,-1,1],Anim1),Anim2), display(op([1,-1,1],Anim1),op([1,-1,1],Anim2),Anim3)], insequence, axes=normal, labels=[x,y,z], view=[0..14,0..14,0..15], orientation=[-45,75]);

                      

Edit.

Download travel.mw

Below is an example of solving your problem for a set of 10 random vectors:

Download GramMat_new1.mw

by the series command.

Example:

N:=10:
series(cos(x)^n, x=0, N);

        1-(1/2)*n*x^2+(-(1/12)*n+(1/8)*n^2)*x^4+(-(1/45)*n+(1/24)*n^2-  (1/48)*n^3)*x^6+(-17/2520*n+7/480*(n^2)-(1/96)*n^3+(1/384)*n^4)*x^8+O(x^10)

on a simple example from a help:

restart;
with(plots): 
A := Matrix([[2, 1, 0, 0, 3], [0, 2, 1, 0, 0], [0, 0, 2, 1, 0], [0, 0, 0, 2, 1], [0, 0, 0, 0, 2]]);
S := [([1, 2, 3, 4, 5]=~StringTools:-Char~(96+~[$1..5]))[]];
sparsematrixplot(A, matrixview, tickmarks = [S, map(t->-lhs(t) = rhs(t) , S)]);

As far as I know, Maple does not implement a step-by-step solution of differential equations. However, you can find out the type of your equation, for example:

DEtools:-odeadvisor(diff(y(x),x)=(x^2+1)*sin(y(x)), y(x));

 Output:                          [_separable]

For details see help on  ?dsolve,details