I made some corrections to your document, and now everything works.

system_solve1.mw

The  createDegree  procedure solves your problem for a general case:

restart;
createDegree:=proc(L::{set,list}, n)
local L1, Degrees, k, Degs, NestedSeq;
L1:=convert(select(t->degree(t)<=n, L), list);
Degrees:=degree~(L1);
k:=nops(Degrees);
Degs:=[seq(floor(n/Degrees[i]), i=1..k)];
NestedSeq:=proc(Expr::uneval, L::list)
local S;
eval(subs(S=seq, foldl(S, Expr, op(L))));
end proc;
[NestedSeq(`if`(`+`(seq(Degrees[j]*(i||j), j=1..k))<=n,`*`(seq(L1[j]^(i||j), j=1..k)), NULL), [seq(i||j=0..Degs[j], j=1..k)])][2..-1];
end proc:

Example of use:
L := {a-2*b, b^2, (a+c)^2}:
for n from 1 to 4 do
createDegree(L,n);
od;

This can be done in several ways. The easiest way is

restart;
plots:-pointplot3d([[1,3,5],[10,30,55],[50,70,25]], style=line, color=red, thickness=2);
 

Example:

sols := [[0,0], [6,0], [6,4], [0,4], [3,6], [6,4], [0,0], [0,4], [6,0]]:
plots:-display(seq(plot(sols[1..n], color=red,thickness=3)$5, n=1..nops(sols)), insequence=true);

Int(Pi*(2*ln(y+sqrt(y^2-4))-2*ln(2))^2, y=2..exp(1)+exp(-1)):
IntegrationTools[Change](%, y+sqrt(y^2-4)=t);  # prompt 1
simplify(value(%));
eval(%, exp(4)-2*exp(2)+1=(exp(2)-1)^2);   # prompt 2
simplify(%) assuming exp(2)-1>0;   # The final result
                         

                                           4*Pi*(exp(1)-4+5*exp(-1))

Addition. Unfortunately Maple itself could not without prompting simplify this:

simplify(sqrt(exp(4)-2*exp(2)+1));   # Should be  exp(2) - 1

Edit.

Here is a simple example with the same error:

restart;
f:=x->diff(g(x), x):
f(0);

     Error, (in f) invalid input: diff received 0, which is not valid for its 2nd argument
 

Maple just substitutes  x=0  and we have  diff(g(0), 0)

A workaround  to avoid the error is usage  D - operator instead of  diff :

restart;
g :=  (x,y)->D[1,1](u1)(x, y)+D[1,2](u2)(x, y);
f :=  y ->g(0, y);
f(0); 
f(1.2);

I understood that the word "simultaneously" means that in every moment  a = b = c :

restart;
eq1:=((diff(f(x),x$3)))+f(x)*diff(f(x),x$2)-a*diff(f(x),x$1)^2=0:
eq2:=(diff(g(x),x$2))+b*f(x)*diff(g(x),x$1)=0:
bc1:=f(0)=0,D(f)(0)=1,D(f)(5)=0,g(0)=0.5,g(5)=0:
sol:=dsolve(subs(a=0.5,b=0.5,{bc1,eq1,eq2}), numeric):
L:=evalf[5]([seq(eval(a*f(x)+b*diff(g(x),x)+c*diff(f(x),x)*g(x),[op(sol(1)),op([a,b,c]=~0.2*i)]), i=0..5)]):
Matrix([[a,b,c,Expr],seq([(0.2*i) $3,L[i+1]], i=0..5)])^%T;

Unfortunately an accurate visualization requires quite a hard work. Here is an example of the visualization of the body of the rotation  a parabola  around the x-axis by disc (washer) method:

f:=x->sqrt(4-x):
with(plots): with(plottools):
A:=plot(f, 0..4, color=black, thickness=2):
B:=plot(f, 0..4, color=yellow, filled=true):
C:=rectangle([2.4,f(2.5)], [2.6,0],color=blue):  # Area element
T:=textplot([[2.5,0,"x",align=below], [2.8,0.7,"f(x)"], [2.5,1.4,dx]], font=[times,roman,14]):
display(C,B,A,T, scaling=constrained, view=[-0.7..4.7,-0.7..2.7]);
A1:=plot3d([f(x)*cos(phi),x,f(x)*sin(phi)], phi=Pi/2..2*Pi, x=0..4, style=surface, color=yellow):
A2:=plot3d([r*cos(phi),0,r*sin(phi)], phi=Pi/2..2*Pi, r=0..2, style=surface, color=yellow):
A3:=plot3d([x,y,0], y=0..4, x=0..f(y), style=surface, color=yellow):
A4:=plot3d([0,y,z], y=0..4, z=0..f(y), style=surface, color=yellow):
B1:=plot3d([f(2.5)*cos(phi),x,f(2.5)*sin(phi)], phi=0..2*Pi, x=2.4..2.6, style=surface, color=blue):
C1:=plot3d([r*cos(phi),2.4,r*sin(phi)], r=0..f(2.5), phi=0..2*Pi, style=surface, color=blue):
C2:=plot3d([r*cos(phi),2.6,r*sin(phi)], r=0..f(2.5), phi=0..2*Pi,style=surface, color=blue):
display(A1,A2,A3,A4,B1,C1,C2, axes=normal, scaling=constrained);  # The body with a cutout for better visibility
Int(Pi*``(f(x))^2, x=0..4)=int(Pi*f(x)^2, x=0..4);  # Calculation of the volume

We have one equation with three unknowns in the real domain. There is an obvious trivial solution (0, 0, 0). Since the equation is homogeneous, if there is a nontrivial solution  (a, b, c) , then  (t*a, t*b, t*c)  also be a solution for any  t>0 . Therefore, the set of all solutions is a cone in R^3 . It is enough to examine the structure of the solutions on the unit sphere. Below we can see from the plot that at the unit sphere we have 3  isolated branches solutions. One solution has been found numerically.

restart;
Eq:=a^2+b^2+c^2+a*b+a*c+b*c-  (a+b-c)*sqrt(2*a*b+a*c+b*c)-(a+c-b)*sqrt(a*b+2*a*c+b*c)-(b+c-a)*sqrt(a*b+a*c+2*b*c):
Eq1 := eval(Eq, [a = cos(phi)*sin(theta), b = sin(phi)*sin(theta), c = cos(theta)]):
plots[implicitplot](Eq1, phi = 0 .. 2*Pi, theta = 0 .. Pi, color = red, grid = [1000, 1000]);
RootFinding[Analytic](eval(Eq1, theta = 0.8), re = 1..2, im = -0.1..0.1);  # a root
evalf(eval(Eq1, [phi = %, theta = 0.8]));  # Check

If I understand your problem  you can use  seq  command.

Example of use:

seq(LinearAlgebra[RandomMatrix](3), i=1..3);

A variation on the theme vv's solution:

restart; 
S := {-x+2*y, -x-y+3*z, 3*x-2*y-z, x+y+z = 1}: 
A, B, C:=map2(eval, [x,y,z], solve~([seq(S minus {S[i]}, i = 1 .. 3)]))[ ]; 
OO := [0, 0, 0]: 
plots:-polygonplot3d([[OO, A, B], [OO, A, C], [OO, B, C]], axes = normal, labels = [x, y, z], labelfont = [times, roman, 16]);

Solid_angle.mw

Edit.

Many integrals (even without parameters ) can not be calculated exactly (symbolically). In this case, the numerical calculation can be used with the specific values of the parameters:

f := GAMMA(phi)*y^(mu*phi-1)*(1-y)^((1-mu)*phi-1)/(GAMMA(mu*phi)*GAMMA((1-mu)*phi)); 
int(eval(ln(1-y)*f, {mu = 1/3, phi = 2}), y = 0 .. 1); 
evalf(%);

v1:=r->rhs(TOV[5])(r)-(rhs(TOV[5])(608529.7)+0.36555);  # The new function
v1(100), v1(1000), v1(100000);  # Examples of use
plot(v1, 0..600000);  # The plot

PDE := {diff(phi(x, t), t) = (diff(phi(x, t), x, x))/S__c-K__r*phi(x, t)+S__r*(diff(theta(x, t), x, x)), diff(u(x, t), t) = diff(u(x, t), x, x)-M^2*(u(x, t)-m*w(x, t))/(m^2+1)-u(x, t)/`&varkappa;`-2*Omega^2*w(x, t)+Gr*theta(x, t)+Gm*phi(x, t), diff(w(x, t), t) = diff(w(x, t), x, x)+M^2*(m*u(x, t)-w(x, t))/(m^2+1)-w(x, t)/`&varkappa;`+2*Omega^2*u(x, t), diff(theta(x, t), t) = lambda*(diff(theta(x, t), x, x))/P__r}:
IBC:={phi(0, t) = 1, phi(9, t) = 0, phi(x, 0) = 0, u(0, t) = t, u(9, t) = 0, u(x, 0) = 0, w(0, t) = 0, w(9, t) = 0, w(x, 0) = 0, theta(0, t) = 1, theta(9, t) = 0, theta(x, 0) = 0}:
P__r:=0.71: lambda:=1.0: S__c:=0.22: K__r:=1.0: S__r:=0.5: m:=0.5: M:=sqrt(10.0): `&varkappa;`:=0.5: Omega:=sqrt(5.0): Gr:=6.0: Gm:=5.0:
sol:=pdsolve(PDE , IBC, numeric, spacestep=0.01);
p1 := sol:-plot(t = 0.3, color = red):
p2 := sol:-plot(t = 0.5, color = gold):
p3 := sol:-plot(t = 0.7, color = yellow):
p4 := sol:-plot(t = 1., color = green):
plots[display]({p1, p2, p3, p4}); 

Edit: spacestep  option has been reduced.

Addition.  The plots of the function phi(x,t) (by default the first function of the system) have been plotted above for the specific values of  t .  If you need the plots of the other functions of the system (for example  u(x,t) ) , you can do it like this:

q1, q2, q3, q4:=seq(eval(u(x,t), sol:-value(t = t0, output = listprocedure)), t0=[0.3, 0.5, 0.7, 1]):
plot([q1,q2,q3,q4], 0..3, color=[red,gold,yellow,green], scaling=constrained, labels=[x, u(x,t)]);

Should be  Md  instead of  Mb :

solve({Fh, Fz, Md, Ms1, Ms2}, {Az, Bz, Dx, Dz, Ma});

      {Az = 44.-1.*Bz-1.*Dz, Bz = Bz, Dx = 1.500000000, Dz = Dz, Ma = 414.-6.*Bz-12.*Dz}

We see that the system has an infinite number of solutions that depend on two parameters  Bz  and  Dz