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)/`ϰ`-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)/`ϰ`+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): `ϰ`:=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

Another way is to use  RealDomain  package:
 

with(RealDomain):
plot(x^(5/3)-5*x^(2/3), x=-3..7, color=red);

In the first example, we use the standard parametrization of the sphere and a rotation. In the second example (the simpler one) the standard parametrization only.

The first example (the upper and lower parts after cutting):

P1:=plot3d([4*cos(phi)*sin(theta),4*sin(phi)*sin(theta),4*cos(theta)], phi=0..2*Pi, theta=0..arccos(-1/2/sqrt(3)), style=surface, color=khaki):
P2:=plot3d([4*cos(phi)*sin(theta),4*sin(phi)*sin(theta),4*cos(theta)], phi=0..2*Pi, theta=arccos(-1/2/sqrt(3))..Pi, style=surface, color=khaki):
P:=plot3d([r*cos(phi),r*sin(phi),-2/sqrt(3)], phi=0..2*Pi, r=0..sqrt(44/3), style=surface, color=khaki):
S1:=plottools[rotate](P1, arccos(1/sqrt(3)), [[0,0,0],[1,-1,0]]):
S2:=plottools[rotate](P2, arccos(1/sqrt(3)), [[0,0,0],[1,-1,0]]):
S:=plottools[rotate](P, arccos(1/sqrt(3)), [[0,0,0],[1,-1,0]]):
plots[display](<plots[display](S1,S) | plots[display](S2,S)>, axes=normal, view=[-5.9..5.9,-5.9..5.9,-4.9..4.9], orientation=[55, 75, 0], scaling=constrained);

The second example (The left and right parts after cutting):

P1:=plot3d([4*cos(theta),4*cos(phi)*sin(theta),4*sin(phi)*sin(theta)], phi=0..2*Pi, theta=0..arccos(5/16), style=surface, color=khaki):
P2:=plot3d([4*cos(theta),4*cos(phi)*sin(theta),4*sin(phi)*sin(theta)], phi=0..2*Pi, theta=arccos(5/16)..Pi, style=surface, color=khaki):
P:=plot3d([5*cos(theta)-2,5*cos(phi)*sin(theta),5*sin(phi)*sin(theta)], phi=0..2*Pi, theta=0..arccos(13/20), style=surface, color=khaki):
plots[display](<plots[display](P1,P) | plots[display](P2,P)>, axes=normal, view=[-5.9..5.9,-5.9..5.9,-4.9..4.9], orientation=[55, 75, 0], scaling=constrained);

All 3 methods are incorrect. Here is the corrected first method:

T := plottools[transform]((x, y, z) ->`if`(z <= 2-y, [x, y, z], [x, 1+y/2-z/2, 1-y/2+z/2])):
plots[display](T(plottools[sphere]([0, 0, 0], 4, style=surface, color=khaki, lightmodel=light4, numpoints = 100000)), scaling = constrained);

Clarification. [x, 1+y/2-z/2, 1-y/2+z/2]  is the projection of the point  [x, y, z]  onto the plane  z = 2 - y

Edit.

func := n -> map(t -> t/content(t), {op(expand((a+b+c)^n))}):

Example of use:
func(3);

                     {a^3, b^3, c^3, a*b^2, a*c^2, a^2*b, a^2*c, b*c^2, b^2*c, a*b*c}

If you have a planar region, use the double integral. See this article

If the boundary of a plane region is specified as a piecewise smooth curve, then to calculate double integrals is advisable to use Green's formula. This technique I've used here to calculate the area of a plane figure.