sort(D^4 + D^3 + D^2 +D + C, D,  descending);

For a solution in  Optimization  package it is useful to use the  initialpoint  option, which can be obtained from the plot:

restart;
plot3d(`if`(y >= -(1/13)*x-85/13 and y <= -(25/14)*x+26/7 and y <= (23/12)*x+89/12, 4*x*y-2*y^2+4*x+y+4, NULL), x=-10..10,y=-10..10, numpoints=40000, axes=normal);
Optimization[Maximize](4*x*y-2*y^2+4*x+y+4,{y >= -(1/13)*x-85/13, y <= -(25/14)*x+26/7, y <= (23/12)*x+89/12}, initialpoint=[x=-5,y=-5]);
Optimization[Minimize](4*x*y-2*y^2+4*x+y+4,{y >= -(1/13)*x-85/13, y <= -(25/14)*x+26/7, y <= (23/12)*x+89/12}, initialpoint=[x=3,y=-5]);

P:=x^8 + 2*x^4 + 1; 
Q:=x^5 + 4*x^3 + x^2 + 3*x +3;
3*P+4*Q mod 5;
P-Q mod 5;
expand(P*Q) mod 5;
R:=rem(P, Q, x, 'q') mod 5;  # remainder
q;  # quotient
is(P=expand(q*Q+R) mod 5);  # Check

Edit.
 

Use a semicolon or a colon after each assignment:

equations:={b[1](2)*a[1](2, 1) = 1/2, b[1](2)*a[1](2, 1)^2 = 1/3, b[1](2)*a[2](2, 1) = 1/2, b[1](2)*a[2](2, 1)^2 = 1/3, b[2](2)*a[1](2, 1) = 1/2, b[2](2)*a[2](2, 1) = 1/2, b[1](2)*(a[3](2, 1)+a[3](2, 2))*a[1](2, 1) = 1/3, b[1](2)*(a[3](2, 1)+a[3](2, 2))*a[2](2, 1) = 1/3, b[1](2)*a[2](2, 1)*a[1](2, 1) = 1/3, b[1](1)*(a[3](1, 1)+a[3](1, 2))+b[1](2)*(a[3](2, 1)+a[3](2, 2)) = 1/2, b[1](1)*(a[3](1, 1)+a[3](1, 2))^2+b[1](2)*(a[3](2, 1)+a[3](2, 2))^2 = 1/3, b[2](1)*(a[3](1, 1)+a[3](1, 2))+b[2](2)*(a[3](2, 1)+a[3](2, 2)) = 1/2, b[1](1)+b[1](2) = 1, b[2](1)+b[2](2) = 1};

variables:={a[1](2, 1), a[2](2, 1), a[3](1, 1), a[3](1, 2), a[3](2, 1), a[3](2, 2), b[1](1), b[1](2), b[2](1), b[2](2), b[3](1), b[3](2)};

solve(equations,variables);

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