Simple procedure  cplot  generalizes  Acer's way and shows also by specific color the points and ranges, in which the function is equal to zero. The colors and thickness of the points, lines and sizes of the points are set by default, but can also be specified by the user.

cplot := proc (f::algebraic, r::range, Colors::list := [blue, red, green], R::realcons := 17, t::`=` := thickness = 2)

local x, L, M, M1, Points, P;

uses plots;

x := indets(f)[1];

L := map(op, {solve({lhs(r) <= x, x <= rhs(r), f}, allsolutions, explicit)});

M := select(s->type(s, `=`) or type(s, `<=`), L);

M1 := select(t->type(t, realcons), map(op, M));

Points := plot([seq([M1[i], 0], i = 1 .. nops(M1))], style = point, symbol = solidcircle, symbolsize = R, color = Colors[3]);

display(plot([piecewise(0 < f, f, undefined), piecewise(f < 0, f, undefined), piecewise(f = 0, f, undefined)], x = r, color = Colors, t), Points);

end proc:

Example:

F := (x+1)*x*(x-2):

f := piecewise(x > -1 and x < 0, 0, F);

cplot(f, -2 .. 3, thickness = 3);

The procedure  cplot  is easily modified to show the intervals of increase, decrease, and constancy of the function.

cplot1 := proc (f::algebraic, r::range, Colors::list := [blue, red, green], R::realcons := 17, t::`=` := thickness = 2)

local x, f1, L, M, M1, Points, P;

uses plots;

x := indets(f)[1];

f1 := diff(f, x);

L := map(op, {solve({lhs(r) <= x, x <= rhs(r), f1}, allsolutions, explicit)});

M := select(s->type(s, `=`) or type(s, `<=`), L);

M1 := select(t->type(t, realcons), map(op, M));

Points := plot([seq([M1[i], eval(f, x = M1[i])], i = 1 .. nops(M1))], style = point, symbol = solidcircle, symbolsize = R, color = Colors[3]);

display(plot([piecewise(0 < f1, f, undefined), piecewise(f1 < 0, f, undefined), piecewise(f1 = 0, f, undefined)], x = r, color = Colors, t), Points);  

end proc:

Example:

F := (x+1)*x*(x-2) - 1:

f := piecewise(x > -1 and x < 0, -1, F);

cplot1(f, -2 .. 3);    # Green color - the points in which the derivative is equal to zero

Another example:

f := 2*sin(t)+3*cos(t)-1:

cplot1(f, -Pi .. 2*Pi, thickness = 3) ;

Edit. Codes have been edited and added another example.

For presentation the different surfaces are painted in different colors:

restart;

A := plots[spacecurve]([x, 0, ln(x)], x = 1 .. exp(1), thickness = 2, color = red):

B := plots[spacecurve]([x, 0, -ln(x)], x = 1 .. exp(1), thickness = 2, color = blue):

C := plots[spacecurve]([exp(1), 0, z], z = -1 .. 1, thickness = 2, color = yellow):

N := 100:

S1 := seq(plot3d([[x*cos(phi), x*sin(phi), ln(x)], [x*cos(phi), x*sin(phi), -ln(x)]], x = 1 .. exp(1), phi = 0 .. 2*Pi*k/N, color = [red, blue], style=surface), k = 1 .. N):

S2 := seq(plot3d([exp(1)*cos(phi), exp(1)*sin(phi), z], phi = 0 .. 2*Pi*k/N, z = -1 .. 1, color = yellow, style=surface), k = 1 .. N):

S := plots[display](A, B, C), seq(plots[display](A, B, C, S1[k], S2[k]), k = 1 .. N):

plots[display](S, insequence=true, scaling = constrained, view = [-4.1 .. 4.1, -4.1 .. 4.1, -1.4 .. 1.4], orientation = [-75, 60], axes = normal);

Your surface is composed of two parts: external and internal. The outer part is indeed a sphere and it covers the inner surface. I made a cut-out and painted in different colors, so that you can see the inner surface. Your body is located between these surfaces.

A := plot3d([x, (-x^2+2*x)*cos(phi), (-x^2+2*x)*sin(phi)], x = 0 .. 2, phi = -Pi .. (1/2)*Pi, color = yellow):

B := plot3d([2*cos(alpha)^2, 2*cos(alpha)*sin(alpha)*cos(phi), 2*cos(alpha)*sin(alpha)*sin(phi)], alpha = 0 .. (1/2)*Pi, phi = -Pi .. (1/2)*Pi, color = cyan):

plots[display](A, B, axes = normal, lightmodel = light4, orientation = [-115, 55]);

Procedures  Y_low, Y_up, Y for each value of  X  are respectively looking  for  Y  for the lower branch, the upper branch and both values. The procedures are based on successive elimination of the variable  theta, beginning with the first equation without  elimination  command. I tried to write everything in one procedure, but problems arose when plotting.

restart;

f := theta->cos(theta)+(2/25)*cos(3*theta):

g := theta->-sin(theta)+(2/25)*sin(3*theta):

Y_low := proc (X)

local m, M, c;

m := minimize(f(theta)); M := maximize(f(theta));

if is(X < m) or is(M < X) then return "No solutions" end if;

c := RealDomain[solve](subs(cos(theta) = y, expand(f(theta) = X)));

simplify(subs({cos(theta) = c, sin(theta) = sqrt(1-c^2)}, expand(g(theta))));

end proc:

Y_up := proc (X)

local m, M, c;

m := minimize(f(theta)); M := maximize(f(theta));

if is(X < m) or is(M < X) then return "No solutions" end if;

c := RealDomain[solve](subs(cos(theta) = y, expand(f(theta) = X)));

simplify(subs({cos(theta) = c, sin(theta) = -sqrt(1-c^2)}, expand(g(theta))));

end proc:

Y:=X->[Y_low(X), Y_up(X)]: 

Examples:

Y_low(1); evalf(%);

Y_up(1); evalf(%);

Y(0);

Y(2);

plots[display](plot(Y_up, -0.5 .. 1, color=red, thickness=2, scaling=constrained), plot(Y_low, -0.5 .. 27/25, color=blue, thickness=2, scaling=constrained));

with(plots):

A := spacecurve([-y^2+4, y, 2-y], y = 0 .. 2, color = red, thickness = 3):

B := plot3d(2-y, x = 0 .. 4, y = 0 .. 2, color = green, style=surface):

C := plot3d([-y^2+4, y, z], y = 0 .. 2, z = 0 .. 2, color = cyan, style=surface):

display(A, B, C, scaling = constrained, axes = normal, view = [-1 .. 4.5, -1 .. 2.4, -1 .. 2.4], lightmodel=light4);

 Addition.  The easiest way of cutting away - to use  filled  option:

restart;

A := plots[spacecurve]([-y^2+4, y, 2-y], y = 0 .. 2, color = red, thickness = 3):

B:=plot3d(2-y, x = 0..4-y^2, y=0..2, color = grey, style=surface, filled, scaling=constrained, lightmodel=light4, axes=normal, view = [-1 .. 4.5, -1 .. 2.4, -1 .. 2.4]):

plots[display](A,B); 

restart;

f :=(theta)->cos(theta)+0.8e-1*cos(3*theta);

g :=(theta)->-sin(theta)+0.8e-1*sin(3*theta);

X^2=f(theta)^2;

Y^2=g(theta)^2;

simplify(expand(%%+%));

solve(%, Y);

You have one nonlinear transcendental equation with 6 unknowns. In general, the set of solutions of similar equations - nonlinear 5-dimensional manifold in the 6-dimensional space. It may turn out to be empty. But in your example, there are infinitely many solutions. The variable  a2  will be regarded as the main variable, other variables will be considered as parameters. Since a2 is included in the equation by linear way, then for each permissible  set of parameters exists not greater than one solution of the equation with respect to a2 matched your criteria. Here is an example of finding one solution (for exact solution all constants should be exact):

restart;

Eq:=50*tan(8*Pi/180)=(a1/(d2-d1))*ln(d2/d1)+(a2/(d3-d2))*ln(d3/d2)+((6/10-a1-a2)/(d4-d3))*ln(d4/d3);

solve(eval(Eq,[d1=5/100, d2=6/100, d3=7/100, d4=9/10, a1=1/10]));

evalf(%);

If you need more solutions, the main problem is the selection of parameters for which and around which there are valid solutions. This can be done in the usual for loop:

N:=0:

for i from 0.01 to 0.06 by 0.01 do

for j from i+0.01 to 1 by 0.1 do

for k from j+0.01 to 1.5 by 0.1 do

for m from k+0.01 to 2 by 0.1 do

for n from 0.1 to 0.5 do

p:=fsolve(eval(Eq,{d1=i,d2=j,d3=k,d4=m,a1=n}));

if p>0 and n+p<0.6 then N:=N+1; L[N]:=[i,j,k,m,n,p] fi;

od: od: od: od: od:

L:=convert(L,list);

Easy modification of  coefff  procedure solves the problem:

coefff:=proc(P, T, t)

local L, H, i, k, s, t1:

L:=[coeffs(P, T, 'h')]: H:=[h]: k:=0:

s:=`if`(t::realcons,t,op(1,t));

if s::realcons and s<>1 then t1:=t/s else t1:=t fi;

for i from 1 to nops(H) do

if H[i]=t1 then k:=L[i] fi:

od:

if t1=t then k else k/s fi;

end proc:

Examples:

f:=a*x^2+b*x*y^3+2;

coefff(f,[x,y],x^2/2);

coefff(f,[x,y],x*y^3);

coefff(f,[x,y],1/10);

The procedure  coefff  returns the coefficient of the monomial  t  in the multivariate polynomial  P (T - the list of variables of the polynomial). The procedure works with both numeric and symbolic  coefficients.

coefff:=proc(P, T, t)

local L, H, i, k:

L:=[coeffs(P, T, 'h')]: H:=[h]: k:=0:

for i from 1 to nops(H) do

if H[i]=t then k:=L[i] fi:

od:

k;

end proc:

Examples of use:

Pol:=3*x^2-a*x*y+1:

coefff(Pol, [x,y], x^2);

coefff(Pol, [x,y], x*y);

coefff(Pol, [x,y], 1);

                 3

                -a

                 1

If I understand the question right, you need to generate all of the matrix of a given size n, which have k units on the diagonal, other elements are 0  (you have n=10). The procedure SM make this.

SM:=proc(n,k)

local P;

uses combinat;

P:=permute([1$k, 0$n-k]);

seq(Matrix(n, p, scan= diagonal), p=P);

end proc:

Example of use:

SM(5,2);

I do not know commands in Maple, which realize an approximation by piecewise constant functions. It is easy to write a procedure that makes it. The  function should be specified in the nodes [x [i], y [i]] . The procedure determines the value of the function at intermediate points as the average of the values in neighboring nodes.

Code of the procedure:

PWCapprox:=proc(L::list)

local n;

n:=nops(L);

piecewise(seq(op([x>L[i,1] and x<L[i+1,1],(L[i,2]+L[i+1,2])/2]), i=1..n-1));

end proc:

Example of use (approximation of sinusoid arch):

n:=20:

L:=[seq([Pi/n*k,sin(Pi/n*k)], k=0..n)]:

A:=plot(L, style=point,color=red):

B:=plot(PWCapprox(L), x=0..Pi, color=red, thickness=2, discont):

plots[display](A, B, scaling=constrained);

See the help in  GraphTheory  package on the commands  DijkstrasAlgorithm , AllPairsDistance, BellmanFordAlgorithm,  ShortestPath

I think that this inequality intended for manual solution. Here is the such solution with partial  use of Maple. Since the function is continuous, then the solution of the inequality reduces to finding the exact roots of the corresponding equation.

A:=18*9^(x^2+2*x)+768*4^((x+3)*(x-1))-5*6^((x+1)^2):

applyop(expand, {seq([i,2,2], i=1..3)}, A);  # We open the brackets in power exponents

algsubs(x^2+2*x=t, %);  # Substitions  t=x^2+2*x

applyop(expand, {[2,2], [3,2]},%);  # Simplification

map(s->s/4^t, %)=0;  # We divide both sides by  4^t

So we come to the quadratic equation  B  for  s=(3/2)^t

B:=18*s^2+12-30*s=0:

solve(B);

                           1, 2/3

It remains to solve two simple equations  (3/2)^(x^2+2*x)=2/3  and  (3/2)^(x^2+2*x)=1  and original inequality:

solve(x^2+2*x=-1), solve(x^2+2*x=0);

solve((x+1)^2*(x+2)*x>=0);

First, we find the rank of a given system of vectors (let rank is n), and then generate all possible sets of n vectors of the original system and the selection of those of them that have a maximum rank.

Example:

restart;

with(LinearAlgebra):

S:=[seq(RandomVector(3, generator=rand(0..3)), n=1..15)];  # The list of 15 random vectors

Matrix(S):

n:=Rank(%);

T:=combinat[choose](S, n):  # The list of all sets of n vectors of S

nops(T);

k:=0:

for t in T do

if Rank(Matrix(t))=n then k:=k+1; B[k]:=convert(t,set) fi;

od:

B:=convert(B, list):  # The list of all bases 

nops(B);  # The number of bases

B[1..10];  # the first 10 bases of B

I do not know in Maple direct commands that solve both problems for multivariate polynomials. Here are some solutions  for the polynomial  P:

P:=a*x^2*y*z+b*x*y*z+c*x^2*y*z+d*x*y^2*z+e*x*y*z;

L:=[coeffs(P, [x,y,z], 't')];

H:=[t];

add(L[i]*H[i], i=1..nops(L));  # Solution of the first problem

for i to nops(H) do  # Solution of the second problem

if H[i]=x^2*y*z then break fi;

od;

L[i];

                                a+c

L[ListTools[Search](x^2*y*z, H)];  # Shorter variant

                                a+c

Addition: If you want to have the members of the polynomial  P  to be sorted, for example, in lexicographical order, then

sort(add(L[i]*H[i], i=1..nops(L)), [x,y,z]);

        (a+c)*x^2*y*z+d*x*y^2*z+(b+e)*x*y*z