@Jjjones98  Yes you are right. Maple returns the value of  A=Sol[1], but for some reason does not return  Sol[2], and so on. But you can still get the values of the rest of the unknowns if you use explicit formulas, for example, Cramer's rule. So the meaning of  B=Sol[2]  will be
Determinant(`<|>`(seq(`if`(k=2,b,A(..,k)),k=1..6)))/Determinant(A);
and so on. Unfortunately, the expressions are still very bulky, and I don’t know how to simplify them and whether this is possible in principle.

Edit.

@Carl Love  In the general case, the result may depend on the order in which the substitutions occur. Compare:

A:=a/x+b*x+1:
S:={a=0, b=2};
eval(eval(A, S), x=0);
eval(A, S union {x=0});

@Yaongyaong  To find the surface area of a paraboloid lying inside a cone, you can use the standard formula through the double integral:

z:=12-x^2-y^2:
int(sqrt(1+diff(z,x)^2+diff(z,y)^2),[y=-sqrt(9-x^2)..sqrt(9-x^2), x=-3..3]);

 Output:                           37*sqrt(37)*Pi*(1/6)-(1/6)*Pi

Another way:

Student:-Calculus1:-SurfaceOfRevolution(12-x^2,x=0..3, axis =vertical, output=value);

Output:                                 (1/6)*(37*sqrt(37)-1)*Pi 

These 2 expressions are not equal:

restart;
A:=Ve*k21*ke*x1+Ve*k21*x1^2+k12*ke^2*x1[t]+2*k12*ke*x1*x1[t]+k12*x1^2*x1[t]+k21*ke^2*x1[t]+2*k21*ke*x1*x1[t]+k21*x1^2*x1[t]+Ve*ke*x1[t]+ke^2*x1[t,t]+2*ke*x1*x1[t,t]+x1^2*x1[t,t]:
B:=x1[t,t]*(x1+k2)^2+(k12+k21)*x1[t]*(x1+ke)^2+Ve*x1[t]*ke+(k21*Ve)*x1*(x1+ke):

expand(A-B);

Output:    -k2^2*x1[t, t]-2*k2*x1*x1[t, t]+ke^2*x1[t, t]+2*ke*x1*x1[t, t]

@acer  We can replace the last 2 lines of your code with one line:

subs(1-isolate(e_nsp2,gamma), isolate(e_nsp1,gamma), set_3);

Maple sequentially performs these 2 substitutions.

@lhasif  Thank you. The surface will look more expressive if you use the following code:

with(Student:-Calculus1):
f := x->(x-1)^2:
g := x->x+1:
A:=SurfaceOfRevolution(f(x), x = 0 .. 3, distancefromaxis = 4, surfaceoptions=[color=khaki, style=patch], output=plot):
B:=SurfaceOfRevolution(g(x), x = 0 .. 3, distancefromaxis = 4, surfaceoptions=[color=khaki, style=patch], output=plot):
plots:-display(A,B, caption=``, lightmodel=light3, scaling=constrained);

For even greater expressiveness, we can make a cutout on this surface in order to clearly see the curve to be rotated and the axis of rotation. But in this case, we must specify the surface and these curves by explicit equations:

f := x->(x-1)^2;
g := x->x+1;
A:=plot3d([[x,(4-f(x))*cos(phi),(4-f(x))*sin(phi)+4],[x,(4-g(x))*cos(phi),(4-g(x))*sin(phi)+4]],x=0..3,phi=-Pi/2..5*Pi/6, style=surface, color=khaki):
B:=plots:-spacecurve([[t,0,f(t)],[t,0,g(t)]], t=0..3, color=red, thickness=3): # The line to be rotated
C:=plots:-spacecurve([t,0,4], t=-1..4, linestyle=3, color=blue, thickness=2): # Axis of rotation
plots:-display(A,B,C, axes=normal, scaling=constrained, view=[-1..4,-4..4,0..8], orientation=[-125,75]);

@Teep You can do it similarly:

Download MaplePrimes_Example_animChart.mw

@Teep    I don’t understand what you call an animated column chart. Write the code for a non-animated chart, I'll try to animate it.

@Carl Love  OP wrote  "...n is natural and k=0..n"

Of course, your interpretation is also possible, but only the OP  can clarify what was meant.

@Carl Love  Of course, it is always useful to have an alternative solution. In my opinion, it’s easier to interpret your line 

B:= eval(L, t= solve(L.v = C.v)); 

(without mentioning the plane) how "the vector  L - C  is perpendicular to the vector  v " and rewrite it as

B:= eval(L, t= solve((L-C).v)); 

@tomleslie  I corrected my answer.

@Zeineb To get a solution, we can take arbitrary vectors, but  n  must be specific. This solution can be easily written as a procedure in which  n  will be one of the parameters. But when using this procedure, we must indicate the specific value of n .

@JTB2222   Similarly.

@ActiveUser 

You wrote a lot of monotonous code, in each line of which you assign a certain expression to  g , and then do some substitution. If you think Maple is calculating something incorrectly, select only this line from the entire code and ask a specific question.

@Zen 

evalindets(g*a*(a*g+g*b+u^2*(Df-Dr))*b/(Df*Dr*(a+b)*u),`+`, t->simplify(t, {a+b=L}));

                             g*a*(L*g+u^2*(Df-Dr))*b/(u*L*Dr*Df)