restart:    
with(IntegrationTools):
M:=Int(lambda(s)*diff(u(s),s,s,s), s);
for n from 1 to 3 do
Parts(M, `if`(n=1, lambda(s), diff(lambda(s), s$(n-1))));  M:=%;
od:
M;

If you specify parameter values, then it is easy to find all the solutions numerically.

Example:

restart;
eq1 := (1/1152)*(81*a1*P^4*b^4*C1^3+162*a1*P^4*b^4*C1*C2^2+648*a1*P^4*b^4*C1*C3^2+432*a1*P^4*b^4*C1*C4^2+864*a1*P^4*b^4*C2*C3*C4+18*a2*P^4*a^2*b^2*C1^3+90*a2*P^4*a^2*b^2*C1*C2^2+90*a2*P^4*a^2*b^2*C1*C3^2+288*a2*P^4*a^2*b^2*C1*C4^2+81*a3*P^4*a^4*C1^3+648*a3*P^4*a^4*C1*C2^2+162*a3*P^4*a^4*C1*C3^2+432*a3*P^4*a^4*C1*C4^2+864*a3*P^4*a^4*C2*C3*C4+288*d1*P^4*b^4*C1+576*d2*P^4*a^2*b^2*C1+288*d3*P^4*a^4*C1+1152*d4*P^4*a^2*b^2*C1+4096*Tc*a^3*b^3*C4)/(b^3*a^3) = 0:

eq2 := (1/1152)*(162*a1*P^4*b^4*C1^2*C2+864*a1*P^4*b^4*C1*C3*C4+81*a1*P^4*b^4*C2^3+432*a1*P^4*b^4*C2*C3^2+648*a1*P^4*b^4*C2*C4^2+90*a2*P^4*a^2*b^2*C1^2*C2+72*a2*P^4*a^2*b^2*C2^3+612*a2*P^4*a^2*b^2*C2*C3^2+360*a2*P^4*a^2*b^2*C2*C4^2+648*a3*P^4*a^4*C1^2*C2+864*a3*P^4*a^4*C1*C3*C4+1296*a3*P^4*a^4*C2^3+432*a3*P^4*a^4*C2*C3^2+2592*a3*P^4*a^4*C2*C4^2+288*d1*P^4*b^4*C2+2304*d2*P^4*a^2*b^2*C2+4608*d3*P^4*a^4*C2+4608*d4*P^4*a^2*b^2*C2-4096*Tc*a^3*b^3*C3)/(b^3*a^3) = 0:

eq3 := (1/1152)*(648*a1*P^4*b^4*C1^2*C3+864*a1*P^4*b^4*C1*C2*C4+432*a1*P^4*b^4*C2^2*C3+1296*a1*P^4*b^4*C3^3+2592*a1*P^4*b^4*C3*C4^2+90*a2*P^4*a^2*b^2*C1^2*C3+612*a2*P^4*a^2*b^2*C2^2*C3+72*a2*P^4*a^2*b^2*C3^3+360*a2*P^4*a^2*b^2*C3*C4^2+162*a3*P^4*a^4*C1^2*C3+864*a3*P^4*a^4*C1*C2*C4+432*a3*P^4*a^4*C2^2*C3+81*a3*P^4*a^4*C3^3+648*a3*P^4*a^4*C3*C4^2+4608*d1*P^4*b^4*C3+2304*d2*P^4*a^2*b^2*C3+288*d3*P^4*a^4*C3+4608*d4*P^4*a^2*b^2*C3-4096*Tc*a^3*b^3*C2)/(b^3*a^3) = 0:

eq4 := (1/144)*(54*a1*P^4*b^4*C1^2*C4+108*a1*P^4*b^4*C1*C2*C3+81*a1*P^4*b^4*C2^2*C4+324*a1*P^4*b^4*C3^2*C4+162*a1*P^4*b^4*C4^3+36*a2*P^4*a^2*b^2*C1^2*C4+45*a2*P^4*a^2*b^2*C2^2*C4+45*a2*P^4*a^2*b^2*C3^2*C4+36*a2*P^4*a^2*b^2*C4^3+54*a3*P^4*a^4*C1^2*C4+108*a3*P^4*a^4*C1*C2*C3+324*a3*P^4*a^4*C2^2*C4+81*a3*P^4*a^4*C3^2*C4+162*a3*P^4*a^4*C4^3+576*d1*P^4*b^4*C4+1152*d2*P^4*a^2*b^2*C4+576*d3*P^4*a^4*C4+2304*d4*P^4*a^2*b^2*C4+512*Tc*a^3*b^3*C1)/(b^3*a^3) = 0:

S:=indets([eq1,eq2,eq3,eq4]) minus {C1,C2,C3,C4};  # The list of the parameters

assign(S=~{$1..nops(S)});  # Assigning numeric values to parameters

S;  # Values of the parameters

evalf([solve({eq1, eq2, eq3, eq4}, {C1, C2, C3, C4})]);  # All the numerical solutions

eval~(5*x, x=~[1,4,10] );

                                    [5, 20, 50]

Of course this is inferior in its capabilities to acer's method, but it is simpler:

Download textplot.mw

solve~(5*x=~[5, 20, 50] );

                               [1, 4, 10]

restart;
m1 := Matrix([[1, 2], [3, 4]]);
m2 := Matrix([[5, 6], [7, 8]]);
'm1'.'m2' = m1.m2;
#  or
'm1'*'m2' = m1.m2;
# or
InertForm:-Parse("m1*m2"):
InertForm:-Typeset(%)=m1.m2;

Edit.

See wiki  https://en.wikipedia.org/wiki/Music_and_mathematics

See:

l := t -> 0.5*tanh(0.5*t):
deq := diff(f(t), t)*l(t)*(diff(f(t), t, t)*l(t)+9.8*sin(f(t)))+diff(l(t), t)*(diff(f(t), t)^2*l(t)-9.8*cos(f(t))+4*(l(t)-0.5)) = 0:
eval(deq, t=0);
eval(%, f(0)=0);

                             -0.500 - 2.450 cos(f(0)) = 0
                                         -2.950 = 0

This is a well-known method. I made a cut for better visibility and chose the same scale on the axes:

f:=r->exp(-(r-3)^2*cos(4*(r-3)));
plot3d([r*cos(phi), r*sin(phi), f(r)], r=1..5, phi=Pi/3..2*Pi, scaling=constrained, axes=normal);

Above is the rotation around y-axis (y-xis pointing upwards).

Addition. Rotation around x-axis (with a cut):

restart;
f:=r->exp(-(r-3)^2*cos(4*(r-3)));
plot3d([r,f(r)*cos(phi),f(r)*sin(phi)], r=1..5, phi=Pi/2..2*Pi, style=surface, scaling=constrained, axes=normal);

Animation the rotation around x-axis:

restart;
f:=r->exp(-(r-3)^2*cos(4*(r-3)));
plots:-animate(plot3d,[[r,f(r)*cos(phi),f(r)*sin(phi)], r=1..5,phi=0..a, scaling=constrained, axes=normal], a=0..2*Pi, frames=90);

rotations.mw

You use the same name  B  to create different objects. List B is not used at all in the code.

Corrected code:

restart;
interface(rtablesize=50);
E:=[E1,E2,E3]:
P:=[0,1,2]:
N:=3:
for b from 1 to 12 by 5 do
F:=(i,j)->
if i=b+P[a] and j=b+P[a] then E[a]
elif i=b+P[a] and j=b+N+2+P[a] then -E[a]
elif i=b+P[a]+N+2 and j=b+P[a] then -E[a]
elif i=b+P[a]+N+2 and j=b+P[a]+N+2 then E[a]
else 0:
end if:
B[b]:=add(i,i=[seq(Matrix(20,F), a=1..3)]);
end do;
H:=[seq(B[b],b=1..12,5)];

Use  seq  command for this:

v := Vector[row]([seq(4*n, n = 1 .. 6)]);

The code will be quite simple if you do not draw a vertical line. It is also easy to depict the curve and the region in different colors:

plots:-display(
   plot(1/ra, x= 0..0.5, 0..10, color=yellow, filled),
   plot(1/ra, x= 0..1, 0..10, color=blue, thickness=2));

I think it will be easier to immediately specify  F  as a function:

F:=(x,y)->f(x)+g(y);
F(x,x);

View -> Show/Hide Contents 

Remove the tick from the "input"

Unfortunately Maple unlike Mathematica weakly solves the problems of an optimization symbolically (exactly). Here are two simple examples:

maximize(x^2+y^2, x=0..1, y=1);  # OK for rectangle region
maximize(x^2+y^2, x=0..1, y=1-x);  # Maple fails for triangular region

This example  z=(1-y^2)/(x^2), x^2 + y^2 <= 1, x>=1/2  is easily solved without Maple by hand. Consider cross-sections of the surface by planes  x=C, where  C  is constants. We have z=(1-y^2)/(C^2)=1/C^2 - 1/C^2*y^2 .  These curves are parabolas (branches down) with vertices at points  (C, 0, 1/C^2) . As  C  is increased, the vertices goes down, so the  maximum is 4  and it reached at the point  x=1/2, y=0 .

Visualization:

plots:-display(seq(plots:-spacecurve([x,y,(1-y^2)/(x^2)],y=-sqrt(1-x^2)..sqrt(1-x^2), color=red, labels=["x","y","z"]),x=0.5..1,0.05));