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

Example with numeric parameters:

with(geom3d):
distance(point(A,1,2,3), plane(p, x+y+z+1=0, [x,y,z]));

                                       7/3*3^(1/2)

Unfortunately it does not work with symbolic parameters:

distance(point(A,a,b,c), plane(p, d*x+e*y+f*z+h=0, [x,y,z])) assuming d^2+e^2+f^2<>0;

       Error, (in assuming) when calling 'geom3d:-plane'. Received: 'unable to define the plane'

Addition. The obvious workaround - just take the general formula and substitute your parameters:

subs({x=a,y=b,z=c}, abs(d*x+e*y+f*z+h)/sqrt(d^2+e^2+f^2));

                      abs(a*d+b*e+c*f+h)/(d^2+e^2+f^2)^(1/2)

If I understand the question correctly, you're trying to plot a partial sums of the corresponding series. The variable  x  can only take discrete values  1, 2, ... . This can be done as follows (on the same plot the sum of the series is constructed, ie 1):

plot([[seq([x, sum(.5^n, n = 1 .. x)], x = 1 .. 10)], 1], x = 1 .. 10, color = [red, blue], style = [point, line], linestyle = 3, symbol = solidcircle, view = [-.5 .. 10, -.1 .. 1.2], size = [800, 300], tickmarks = [10, default]);

 Edit. The code was edited.

restart;

Sys:=[x1*a+x2*b+x3*c+x4*d=0,

y1*a+y2*b+y3*c+y4*d=0,

z1*a+z2*b+z3*c+z4*d=0,

t1*a+t2*b+t3*c+t4*d=0];

Var:=[a,b,c,d];

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

A . Vector(Var) = v;  # Your system in matrix form

Addition: If you want to prevent the immediate calculation of the left (for example, for presentation) you can write as follows. On subsequent calls the calculation will be done:

subs(X=A, X.Vector(Var) = v);

%;

1) Your task is equivalent to the problem of extraction coefficients of a multivariate polynomial. First, you generate the list of all key variables (the list [C] as Rouben wrote).  But no need to write  u^k, just write  u  in the Rouben's code for  C.

2) If you need to get the coefficient of the specific monomial (monomial is the product of the key variables in some degrees), the following simple procedure can help you. Procedure  coefff returns the coefficient of the monomial  t  in the multivariate polynomial  P  (parameter  T  is the list of the key variables):

coefff:=proc(P,T,t)

local L,H,i,k:

L:=[coeffs(P,T,'h')]: H:=[h]: k:=0:

for i from 1 to nops(H) do

if H[i]=t then k:=L[i] fi:

od:

k;

end proc:

Example of use:

coefff(a*u[x]^2+b*u[x]*u[y]+c*u*u[y]^2, [u,u[x],u[y]], u[x]*u[y]);
                                                     b

Two methods of the example :

A:=<x,x^2; x^3,x^4>; B:=<cos(x),-sin(x); sin(x),cos(x)>;

C:=(A^+).B;

C1:=map(diff,C,x);  # Direct elementwise differentiation

C2:=map(diff,A^+,x).B + (A^+).map(diff,B,x);  # Differentiation by OP's formula

A^+  is a transposed matrix to  A .

Addition:  In recent versions Maple, you can use a more compact syntax:

C1:=diff~(C,x);