Use  seq  command instead of for loop:

L:=[-3, -2.5, 0, 1, 3/2, Pi, sqrt(5)];

seq([s, type(s,nonnegint)], s=L);

A:=Matrix([[1,2,3], [4,5,6]]);

A[..,1..2];

I have changed the last two lines of Preben's code:

1) With the command  simplify  with option  zero , I removed all the zero imaginary parts.
2) The vector of eigenvalues was converted to the list.
3) The matrix of eigenvectors converted into a sequence of explicit column-vectors

convert(simplify(L, zero), list);   #The list of eigenvalues

seq(simplify(V, zero) [ .. , i], i = 1 .. 12);   #The  eigenvectors corresponding to the eigenvalues, in order

by only the single  plot  command:

plot([x^2-4, [[-2, 0],[2, 0]]], x = -2.9 .. 2.9, style = [line, point], color = [blue, red], thickness = 2, symbol= solidcircle, symbolsize = 20, scaling = constrained);

with(geometry):

Eq:=x=9*t^2 + 3*t  + 367,y=3*t^2 + 2*t + 122:

eliminate({Eq}, t);

conic(c,op(%[2]),[x,y]);

detail(c);

plot([rhs~([Eq])[],t=-0.5..0.1], color=red,thickness=3);

Plotting on the basis of parametric equations (without additional options) provides higher quality graphics than with the implicit equation.

Of course, this method is more laborious than the above, but in some situations it may be beneficial:

applyop(collect, [1,2,1,3,1], phi, cos(2*theta));

Simple procedure  WideMatrix  due to the introduction of empty columns moves apart the columns so that they are clearly visible.

WideMatrix := proc(A::Matrix)

local m, n;

m:=op([1,1],A); n:=op([1,2],A);

Matrix(m,2*n-1, (i,j)->`if`(j::odd,A[i,(j+1)/2],``));

end proc:

Example of use:

A := LinearAlgebra[RandomMatrix](5);

WideMatrix(A);

To construct a surface of revolution, you can use a parametric form of a surface.

The surface of the revolution around the polar axis (the axis Ox):

p := 1-cos(theta):

plot3d([p*cos(theta), p*sin(theta)*cos(phi), p*sin(theta)*sin(phi)], theta = 0 .. 2*Pi, phi = 0 .. 2*Pi, axes = framed, labels = [x, y, z], scaling = constrained);

The surface of the revolution around the axis Oy:

plot3d([p*cos(theta)*cos(phi), p*sin(theta), p*cos(theta)*sin(phi)], theta = 0 .. 2*Pi, phi = -(1/2)*Pi .. Pi, axes = framed, labels = [x, y, z], scaling = constrained);

The cut made to make it easier to see the surface. 

It seems you want it:

A:=Matrix(1,2, {(1,1)=Matrix([[a,b],[c,d]]), (1,2)= Matrix([[e,f],[g,h]]) });

B:=Matrix([A[1,1], A[1,2]]);

Above the polar equation is obtained for the ellipse centered at the origin. In connection with applications in physics interesting polar equation of an ellipse whose focus is at the origin (for definiteness - left focus). For details see https://en.wikipedia.org/wiki/Ellipse .  That's the conclusion of such an equation in Maple:

subs({x = r*cos(phi), y = r*sin(phi)}, (x-sqrt(a^2-b^2))^2/a^2+y^2/b^2 = 1);

solve(%, r);

simplify([%]) assuming a > 0, b > 0;

simplify(subs(b^2 = a^2-f^2, %)) assuming f > 0;

Eq := r = simplify(subs(f = e*a, %[1]));  # Polar equation with the left focus in the origin 

plots[polarplot](eval(rhs(Eq), {a = 3, e = 2/3}), phi = 0 .. 2*Pi, color = red, thickness = 3, scaling = constrained);

In Maple, you can operate with symbolic matrices specified sizes.

Example:

A:=Matrix(3, symbol=a);

AT:=A^+;  # The matrix transpose of the matrix A

A.AT;

Use the following compact syntax to input the system:

Var:=[a,b,c,d,e,f,g,h]:  # List of unknowns

C:=[x,y,z,t,s,p,q,u]:   # List of coefficients

Sys:=[seq(add(cat(C[i],j)*Var[j], j=1..8), i=1..8)];

(A, v) := LinearAlgebra[GenerateMatrix]( Sys, Var );

A . Vector(Var) = v; 

seq(op([a[i],b[i]]), i=1..3);

       a[1], b[1], a[2], b[2], a[3], b[3]

Another way:

f := i -> (a[i], b[i]):
seq(f(i), i=1..3);

Construct a straight line in Maple can be a variety of ways. Сarl showed one way in the package  geom3d. To this end, I usually use the package  plots. In this package, there are 2 ways to build a  line in 3D. The first way is by spacecurve  command if a line defined by parametrically  x=x0+a*t, y=y0+b*t, z=z0+c*t . The second way if a line defined as the intersection of two planes. 

In the example below, the same line is constructed by these methods:

with(plots):

spacecurve([2-t,2-2*t,t], t=0..1, color=red, thickness=3, axes=normal, orientation=[25,60], labels=[x,y,z], view=[0..2,0..2,0..1], scaling=constrained);

intersectplot(x-y-z=0, x+z-2=0, x=0..2,y=0..2,z=0..1,color=red, thickness=3, axes=normal, orientation=[25,60], labels=[x,y,z], scaling=constrained);

The first method gives a higher quality of plotting.

At first I clearly decided your system for the derivatives. We see that it is divided into two systems. The solution of the first system is plotted. To curve was better seen, I have increased the range in all variables. Since your system is not autonomous, there are no any arrows.

ode1:=a(t)*(diff(a(t), t))+c(t)*(diff(c(t), t))=3*t:

ode2:=a(t)*(diff(a(t), t))+a(t)*(diff(a(t), t))*b(t)*(diff(b(t), t))*c(t)*(diff(c(t), t))=3*t:

ode3:=a(t)*(diff(a(t), t))+a(t)*(diff(a(t), t))*c(t)*(diff(c(t), t))+a(t)*(diff(a(t), t))*b(t)*(diff(b(t), t))*c(t)*(diff(c(t), t))=2*t:

subs({diff(a(t), t)=A,diff(b(t), t)=B,diff(c(t), t)=C}, {ode1,ode2,ode3}):

solve(%,{A,B,C});

allvalues(%);

Sys:=subs({A=diff(a(t), t),B=diff(b(t), t),C=diff(c(t), t)},%[1]);

DEtools[DEplot3d](Sys,[a(t),b(t),c(t)],t=0..3,a=0..3,b=0..3,c=0..3,[[a(1)=1,b(1)=2,c(1)=3]],scene=[a(t),b(t),c(t)], linecolor=red, axes=normal, orientation=[30,65]);