@Daniel Maxin  You can use  listcontplot  command, but for this you have to rewrite your data so that it was a two-dimensional list consisting only of the  z-values corresponding to the points  (x[i], y[j]), taking by "layers". There is another problem. Maple draws the contour curves with the marks on the coordinate axes are not real values  (x[i], y[j]), and their numbers  i=1..m, j=1..n . This problem can be easily solved using a linear transformation which is the compression to the point  (1, 1), then the parallel translation.

For example, I took the values  (x[i], y[j])   in the rectangle  x=1..4  and  y=2..4  with step  0.2. Third coordinate z for these points I took by formula  z=x^2+y^2 . So I got the list  Data . Then  plots[listcontplot]  command for  Data . The color intensity increases from pink to red with increasing  z-values. Default Maple draws 8 contour lines. Their numbers and values you can adjust by levels  option.

X:=[seq(x, x=1..4, 0.2)];

Y:=[seq(y, y=2..4, 0.2)]; 

Data:=[seq([seq(x^2+y^2, x=X)], y=Y)];

P:=plots[listcontplot](Data, filledregions=true, coloring=[pink,red]):

P1:=plottools[translate](plottools[homothety](P, 1/5, [1,1]), 0, 1):

plots[display](P1, view=[-0.5..5,-0.5..5]); 

In older versions of Maple in which there is no  parametric  option, it is useful to factorize the polynomial , then the solution is reduced to the study of the discriminant of quadratic trinomial:

factor(b^4-(2-d)*b^2-2*d*b+1+d);

discrim(op(2, %), b);

Because  for  d>0  the discriminant  of  b^2+2*b+d+1  is negative, then  b^4-(2-d)*b^2-2*d*b+1+d>0  for all  b<>1 .  So should be solved in this way, if you need a manual solution.

For automated solution you can simply write

solve(b^4-(2-d)*b^2-2*d*b+1+d > 0, b) assuming d > 0, d <= 1;

It's probably incorrect syntax. See the example:

plott(x^2, x=0..1);

If my example does not help, please provide your code for more accurate response.

Addition. If at first sight the command looks right, check that all the letters typed in English layout.

Probably this is an anomaly of specific version of Maple. For example I have no problem in classic worksheet Maple 12:

plot(1/cos(t), t=0..Pi/3, thickness=2, coords=polar);

Addition. I remember that the problems with polar plots were in Maple 16.00

1. No need to call  plots:-shadebetween  command, as if the lower plane  z = 0 , older option  filled=true  does the same thing.

2. For an accurate representation of the shape of the body  scaling=constrained option should be used.

3. To improve the appearance  user lighting was used (by right clicking).

plot3d(x^2-y^2+4, x = -1 .. 1, y = 0 .. 2, filled = true, axes = normal, view = [-1.4 .. 1.4, 0 .. 2.4, 0 .. 5.4], scaling = constrained, orientation = [50, 70]);

Eq := 544*X-228412+1836*Y-1296*Z = 0:

d := igcd(coeffs(lhs(Eq)));

d*``(lhs(Eq)/d) = 0;

plots[implicitplot](1=x^2+(y-1)^2, x=-2..2, y=-1..3, color=red, thickness=2, view=[-1.5..1.5,-0.5..2.5], gridrefine=3);

For plotting you must specify   h  and  b .

Example:

h:=-x:  b:=(x+3)^3:

f:=piecewise(-3 < x and x< -2, h, -2 < x and x< -1, b):

plot(f, x=-3..-1, discont, thickness=2, view=[-3.5..0.5,-0.5..3], scaling=constrained);

We need to find a point on the curve, in which  the normal passes through original point. The following code solves the problem for the original point  M=[0,-2] :

restart;

f:=x->1/x:

Eq:=(y=f(x0)-1/D(f)(x0)*(x-x0)):

M:=[0,-2]:

x0:=fsolve(eval(Eq, [x,y]=~M), x0=0..infinity);

A:=plot([1/x, f(x0)+D(f)(x0)*(x-x0), [t,eval(rhs(Eq),x=t), t=0..x0]], x=-2..4, -3..3, color=[red,blue,green], thickness=[2,1,2], discont):

B:=seq(plottools[disk](p,0.03,color=yellow), p=[M,[x0,f(x0)]]):

C:=plots[textplot]([[0.2,-2,"M"],[x0+0.2,f(x0)+0.1,P]], font=[TIMES,ROMAN,18]):

plots[display](A,B,C);  # Visualization - the point P is the projection of the point M onto the curve  y=1/x, x>0

 Addition: In general, the solution is not unique. But in this example, for each of the branches of the hyperbola (due to the convexity) solution is unique.

Acer method probably the shortest. Here are more 3 ways:

X := [1, 3, 8, 9, 10]:

`+`(seq(X[i]*ln(X[i]), i=1..nops(X)));

# Or

`+`(op(X*~(ln~(X)))); 

# Or in for loop

S := 0: for x in X do S := S+x*ln(x) end do: S;

                         21*ln(3)+24*ln(2)+10*ln(10)

                         21*ln(3)+24*ln(2)+10*ln(10)

                         21*ln(3)+24*ln(2)+10*ln(10)

To see the point of intersection, you can simply reduce the ranges of the axes at standard numpoints:

restart;

convert([92*x/790+18*x*(1-y)/1000=0.125,ln(46.59/x)=1/1.5*(ln(0.553/y)+2.5*ln((1-y)/(1-0.553)))], rational);

fsolve(%);

plots:-implicitplot(%%, x=0..1, y=0..0.05, thickness=2, color=[blue,red]);

Your system can easily be solved exactly. In general case, it will have  2^3 = 8  solutions. Since the orthogonal matrix  RT has determinant  1, the matrix  mat  should also be orthogonal with determinant  1. For solution we need to select 3 independent elements of the orthogonal matrix  mat (I took  mat[1,2], mat[1,3], mat[2,3]). Here are all the solutions for the specific orthogonal matrix  mat  with full check of one solution:

restart;

Rx := Matrix(3,3, [1,0,0,0,cos(a),-sin(a),0,sin(a),cos(a)]):   

Ry := Matrix(3,3, [cos(b), 0, -sin(b), 0, 1, 0, sin(b), 0, cos(b)]):  # Corrected matrix

Rz := Matrix(3,3, [cos(c), -sin(c), 0, sin(c), cos(c), 0, 0, 0, 1]):

RT := Rx.Rz.Ry:

mat:=<1/3,2/3,-2/3; -2/3,2/3,1/3; 2/3,1/3,2/3>;  # The specific orthogonal matrix with det=1

assign(seq(seq(g[i,j]=(RT[i,j]=mat[i,j]), j=1..3), i=1..3));

Sol:=allvalues([solve({g[1,2], g[1,3], g[2,3]})]);  # Found solutions

eval([seq(seq(g[i,j], j=1..3), i=1..3)], Sol[1,1]);  # The check of  Sol[1,1]  for all equations of the system

Additions.  1) You can verify that  Sol[4,1]  is also a solution. The rest of the solutions are not solutions of the whole system.

2) The code was edited.

Your procedure works for me, but the same can be written a much shorter.

Your way:

dList:=proc(a::list, b::list)

    local res::list, i::integer, n::integer;

    description "add list elements and diff";

    n:=min(nops(a), nops(b)) ;

res:=[seq(0, i=1..n)];

for i from 1 by 1 to n do

     res[i]:=diff(a[i]+b[i],x)

end do;

res;

end proc:

list1:= [x^3 , 2*x , 3]:

list2:= [x^5, 2*x^4 , 2*x]:

dList(list1,list2);

                    [3*x^2+5*x^4, 2+8*x^3, 2]

The shorter method (without user's procedure):

list1:= [x^3 , 2*x , 3]:

list2:= [x^5, 2*x^4 , 2*x]:

map(diff, list1+list2, x);

                     [3*x^2+5*x^4, 2+8*x^3, 2]

 Maybe your error is caused by previous assignments. Make restart before your procedure.

Everything to be beautiful, I changed all the colors and ranges for the axes. Instead of  Student package I just wrote the simple equations for the tangent and normal:

f := x-> x^2-2*x:

x0 := 1.5:

T := (D(f))(x0)*(x-x0)+f(x0):  # The tangant

N := -(x-x0)/(D(f))(x0)+f(x0):  # The normal

P := plottools[point]([x0, f(x0)], symbol = circle, symbolsize = 15):  # The point

Curves := plot([f(x), T, N], x = -4 .. 5, -4 .. 5, color = [red, blue, green], thickness = 2, legend = [curve, tangent, normal]):

plots[display](Curves, P, scaling = constrained);

You will have more opportunities to vary the shape of the knot, if you use it's parametric equations and  plots[tubeplot]  command. You are setting the parameters  rho>1  and coprime integers  p  and  q .

Example:

restart;

r := cos(q*phi)+rho:  x := r*cos(p*phi):  y := r*sin(p*phi):  z := -sin(q*phi):  p := 3:  q := 7:  rho := 3:

plots[tubeplot]([x, y, z], phi = 0 .. 2*Pi, radius = 0.3, style = surface, numpoints = 3000, scaling = constrained, axes = none);