If you want to build several curves defined by explicit equations in the same plot, use square brackets (a list of curves). Here is an example of plotting all 4 basic trigonometric functions:

restart;
plot([sin(x), cos(x), tan(x), cot(x)], x=-Pi/2..5*Pi/2, y=-4..4, color=[red, blue, green, gold], size=[800,500], scaling=constrained,  discont);
 

You can easily do this one line without repeating:

restart;
f(x):=exp(-x+3.):
T__4:=unapply(convert(taylor(f(x), x=7, 5), polynom), x);  # T__4 as an operator
T__4(3);

or (without unapply):

restart;
f(x):=exp(-x+3.):
T__4:=convert(taylor(f(x), x=7, 5), polynom);  # T__4 as an expression
eval(T__4, x=3);

To calculate the values of the derivative at individual points use  unapply  command or  D  operator:

KK := unapply(diff(x^2, x), x);
KK(2);  # or
D(x-> x^2)(2);

We just write the variable  s  as an index:

sys :=
  diff(x[s](t),t) = y[s](t),
  diff(y[s](t),t) + x[s](t) = 0;

ic := x[s](0) = a(s),
      y[s](0) = b(s);
dsolve({sys, ic});

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