1) Go into the site mapleprimes.com.

2) Log in.

3) In the upper right corner there is the word Account. Click on it.

4) ...

Painting individual faces can be set by  plottools[polygon]  command. With regard to edges by  plottools[line]  command.

Example:

with(plottools):

cube := cuboid([0, 0, 0], [1, 1, 1], color = black, style = wireframe):

face1 := polygon([[0, 0, 0], [1, 0, 0], [1, 0, 1], [0, 0, 1]], color = yellow):

face2 := polygon([[0, 0, 0], [0, 0, 1], [0, 1, 1], [0, 1, 0]], color = blue):

face3 := polygon([[0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 0, 0]], color = green):

plots[display](cube, face1, face2, face3, scaling = constrained, axes = boxed, thickness = 3, orientation = [40, 75]);

Your question needs to be clarified. You can plot only a function of one or two variables. Therefore:

1) The constants  C1, C3, C4, n  should be specified.

2) You wrote that the numbers  b[m], m=1..n   are the roots of the equation  C3*sin (beta/2) =beta*cos (beta/2). But this equation has an infinite number of zeros. I assumed that you are interested in the first  n  positive roots. This equation can be rewritten as  beta=C3*tan(beta/2)  and it becomes clear that if  0<C3<2 then this equation in each interval  0..2*Pi, 2*Pi..4*Pi, 4*Pi..6*Pi  and so on has only one root.

Example:

restart;

W:=(tau, eta)->tau*exp(C4)/C1-exp(C4)*(exp(C3-C3*eta)/C3^2-exp(C3)/C3^2+tau/C1+eta/C3)+exp(C4-C3*eta/2)*sum(16*beta[m]*sin(beta[m]*eta/2)*exp(-(beta[m]*beta[m]+C3^2)*tau/(4*C1))/((beta[m]*beta[m]+C3^2)*( beta[m]*beta[m]-C3^2)^2),m=1..n):

C1, C3, C4, n:=1, 1, 1, 10:

beta:=convert(`minus`({RootFinding[Analytic](sin((1/2)*beta)-beta*cos((1/2)*beta), re = 0 .. 2*n*Pi, im = -1 .. 1)}, {0.+0.*I}), list);

unapply(W(tau, eta), tau, eta);

plot3d(W(tau, eta), tau=0..4, eta=0..4, axes=normal);

You can also construct graphs of functions of one variable by fixing the value of another variable, such as

plot(W(1,eta), eta=0..4);

or

plot(W(tau,2), tau=0..infinity);

Use  is  command:

is(1.551691464 < (5/6)*Pi);

             true

See  FunctionAdvisor  command.

Example:

L:=[fsolve(x^8-4*x^7+6*x^6-4*x^5+x^4+x^2+2*x+1, x, complex)];

seq(plottools[disk]([Re(L[i]), Im(L[i])], 0.03, color=red), i=1..nops(L)):

plots[display](%);

2 vectors are collinear if and only if the matrix composed of their coordinates has rank<2. This property is an equivalence relation, and therefore the whole set can be divided into classes according to this property. The procedure  P  produces a partition of S and of each class takes one  first element.

P:=proc(S)

local f, L, M;

L:=convert(S, list);

f:=(x,y)->`if`(LinearAlgebra[Rank](Matrix([x,y]))<2, true, false);

M:=[ListTools[Categorize](f, L)];

print(M);

{seq(M[i,1], i=1..nops(M))};

end proc:

Your example:

S:={[-2, 1, 1, -2, 1, 1], [-2, 1, 1, 2, -1, -1], [0, -1, 1, 0, 1, -1], [0, 1, -1, 0, -1, 1]}:

P(S);

[[[-2, 1, 1, -2, 1, 1]], [[-2, 1, 1, 2, -1, -1]], [[0, -1, 1, 0, 1, -1], [0, 1, -1, 0, -1, 1]]]
            {[-2, 1, 1, -2, 1, 1], [-2, 1, 1, 2, -1, -1], [0, -1, 1, 0, 1, -1]}

If the system contains parameters, the  LinearSolve  command does not give the correct result, even if all the dimensions are specified.

Example:

A:=Matrix([[a,1], [1,b]]):

B:=Vector([1,2]):

LinearAlgebra[LinearSolve](A, B);

For example, if  a=1/2  and  b=2, the result is incorrect.

with(plots):

D124 := implicitplot3d(max(-x+y+z, x-y+z, x+y-z) = 1, x = 0 .. 1, y = 0 .. 1, z = 0 .. 1, style=surface, numpoints = 50000):

D134 := implicitplot3d(max(-x-y-z, -x+y+z, x-y+z) = 1, x = -1 .. 0, y = -1 .. 0, z = 0 .. 1, style=surface, numpoints = 50000):

D123 := implicitplot3d(max(-x-y-z, -x+y+z, x+y-z) = 1, x = -1 .. 0, y = 0 .. 1, z = -1 .. 0, style=surface, numpoints = 50000):

D234 := implicitplot3d(max(-x-y-z, x-y+z, x+y-z) = 1, x = 0 .. 1, y = -1 .. 0, z = -1 .. 0, style=surface, numpoints = 50000):

display(D124, D134, D123, D234, axes=normal, view=[-1.7..1.9, -1.7..1.9, -1.7..1.9], orientation=[30, 50], lightmodel=light4);

When you write  contours=n  then maple  finds the corresponding values ​​of the function by the formulas

seq(m+(1/(n+1))*(M-m)*k,  k = 1 .. n) , 

in which  m  and  M  are the maximum and minimum values of the  function in the specified region.

Your example:

with(plots):

z := x^5-y^3:

m := minimize(z, x = -5 .. 5, y = -4 .. 3); M := maximize(z, x = -5 .. 5, y = -4 .. 3);

L := [seq(evalf[4](m+(1/11)*(M-m)*k), k = 1 .. 10)];

A := textplot([seq([surd(L[k]+(-3+.5*k)^3, 5), -3+.5*k, L[k], align = right], k = 1 .. 10)]):

B := contourplot(z, x = -5 .. 5, y = -4 .. 3, axes = boxed, filledregions = true, coloring = [yellow, red], contours = 10):

display(A, B);

A nonlinear equation can have any number of zeros from 0 to infinity. The  evalf  command finds only one root
Therefore,  firstly it is  useful to examine the properties of the function  s->tanh(s)-C3*s. Here are
two basic properties:
 
1) It is odd.
2) On the positive s-axis the function is strictly convex upward.

From these properties it follows that if  0 <C3 <1, then there exists a unique positive root  s> 0. Altogether there will be three roots  -s, 0, s .

The solution of this equation can be issued in the form of a procedure in which the formal arguments are  C3  and  N (number of digits in the result).

Sol := proc (C3, N)

 local sol, h, A, B;

Digits := N;

if 1 <= C3 or C3 <= 0 then print(`One solution`);  print(s = 0);

plot([tanh(s), C3*s], s = -1.5 .. 1.5, -1.5 .. 1.5, color = [red, blue]) else sol := fsolve(tanh(s) = C3*s, s = 0 .. infinity);

print(`Three solutions`); print(s[1] = -sol, s[2] = 0, s[3] = sol);

h := eval(tanh(s), s = sol);

A := plots:-textplot([[-sol, 0, s[1]], [0, 0, s[2]], [sol, 0, s[3]]], align = [left, above]);

B := plot([tanh(s), C3*s, [-sol, t, t = -h .. 0], [sol, t, t = 0 .. h]], s = -sol-1 .. sol+1, thickness = [2, 2, 1, 1], color = [red, blue, black, black], linestyle = [1, 1, 2, 2], `if`(sol < 5, scaling = constrained, NULL));

plots:-display(A, B);

end if;

end proc;

Example:

Sol(0.18, 20);

f := x^5+a*x^4+b*x^3+c*x^2+d*x+e;
algsubs(x^3=1, f);

                   f := x^5+a*x^4+b*x^3+c*x^2+d*x+e
                            (1+c)*x^2+(a+d)*x+b+e

Procedure called  P  is plotting  your function so that x-axis to be plotted  at y=y0 . Parameters  a  and  b  specify the range on x-axis.

P := proc (expr, a, b, y0)

local m, M;

m := min(minimize(expr-y0, x = a .. b), -1);

M := max(maximize(expr-y0, x = a .. b), 1);

plot(expr-y0, x = a .. b, y = m .. M, tickmarks = [default, [seq(i = i+y0, i = floor(m) .. ceil(M))]]);

end proc:

Example:

P(x^2-4, -2, 2, -5);

Edited. The x-axis can be anywhere: below, above or to cross the plot of the function. 

I copied your expression without the last term, as it is trimmed in your message.

A:=-4*cos((1/5)*Pi)*m^3-4*cos((1/5)*Pi)*n^3-4*cos((1/5)*Pi)*p^3-4*cos((1/5)*Pi)*q^3+4*cos((2/5)*Pi)*m^3+4*cos((2/5)*Pi)*n^3+4*cos((2/5)*Pi)*p^3+4*cos((2/5)*Pi)*q^3-10*cos((1/5)*Pi)*p*q^2+10*cos((2/5)*Pi)*m^2*n+12*cos((2/5)*Pi)*m^2*p+12*cos((2/5)*Pi)*m^2*q+12*cos((2/5)*Pi)*m*n^2+10*cos((2/5)*Pi)*m*p^2+12*cos((2/5)*Pi)*m*q^2+12*cos((2/5)*Pi)*n^2*p+10*cos((2/5)*Pi)*n^2*q+12*cos((2/5)*Pi)*n*p^2+12*cos((2/5)*Pi)*n*q^2+12*cos((2/5)*Pi)*p^2*q+10*cos((2/5)*Pi)*p*q^2-10*cos((1/5)*Pi)*m^2*n-12*cos((1/5)*Pi)*m^2*p-12*cos((1/5)*Pi)*m^2*q-12*cos((1/5)*Pi)*m*n^2-10*cos((1/5)*Pi)*m*p^2-12*cos((1/5)*Pi)*m*q^2-12*cos((1/5)*Pi)*n^2*p-10*cos((1/5)*Pi)*n^2*q-12*cos((1/5)*Pi)*n*p^2-12*cos((1/5)*Pi)*n*q^2-12*cos((1/5)*Pi)*p^2*q-24*cos((1/5)*Pi)*m*n*p-24*cos((1/5)*Pi)*m*n*q-24*cos((1/5)*Pi)*m*p*q-24*cos((1/5)*Pi)*n*p*q+24*cos((2/5)*Pi)*m*n*p+24*cos((2/5)*Pi)*m*n*q+24*cos((2/5)*Pi)*m*p*q:

simplify(convert(A, radical));

Example:

M := plottools:-arc([0, 0], 1, (1/6)*Pi .. 3*Pi*(1/4), color = blue, thickness = 4):

plots:-display(plot(op(M)[1], color = pink, filled, scaling = constrained), M);