@Carl Love  This probably only works in recent versions of Maple.

In Maple 2018.2 :

[for i from 11 to 20 do i^2 od];

     Error, reserved word `for` unexpected

@nm   Try  color="DarkSlateGray"  instead of color=gray .

You should present all of this in text form using Maple syntax. Why should someone do this for you?

@Teep  Of course, all this is easy to do on one plane:

restart; with(plots): 
z := (x, y)->(x^2+y*x+y^2)*sqrt((1-y^2)^2/(x+y)^2-(x-y)^2): 
Colors := [yellow, orange, red, green, blue, khaki, violet, cyan,pink,magenta]: 
plot([seq([x, z(x, y), x = 0 .. 1], y = 0 .. 0.9, 0.1)], color = Colors, thickness = 2, legend = [seq(y = C, C = 0 .. 0.9, 0.1)], labels = [x, z]);

Edit.    

@Carl Love  Thank you. I do not see much need to use polar axes, because a large number of additional lines and inscriptions only interfere with better viewing the curve itself.

@JAMET It will be a regular segment:

plot([r*cos(Pi/4),r*sin(Pi/4),r=1..2], color="Green", thickness=3, view=[0..1.5,0..1.5]);

@emendes  Then I do not understand you. You wrote "...I just need subsets with a chosen number of elements..". I showed how to generate all subsets of a given size. You write that this is not what you need. Can you more accurately formulate your goal?
Here's an example of another way to solve the problem. We generate all triples of numbers in the range 1..10, the sum of which does not exceed 10:

seq(seq(seq(`if`(i+j+k<=10,{i,j,k},NULL), k=j+1..10), j=i+1..9), i=1..8);

   {1, 2, 3}, {1, 2, 4}, {1, 2, 5}, {1, 2, 6}, {1, 2, 7}, {1, 3, 4}, {1, 3, 5}, {1, 3, 6}, {1, 4, 5}, {2, 3, 4}, {2, 3, 5}

See also the commands in the  Iterator  package for your purposes.

@SIMBA  It is probably impossible to obtain the indefinite integral 

int(tanh(x)/sqrt(x^2+1), x);

 in closed form.

@Al86  In order for the equality  f'(x1)=f'(x2)=(f(x1)-f(x2))/(x1-x2)  to be fulfilled, you need to find such 2 points on the graph so that the line passing through these 2 points touches the graph at these points. The solution is below.

restart;
f:=x->a/2*x*(1-x)+(x/4)*ln(x)+(1-x)*ln(1-x):
a:=3:

fsolve({D(f)(x1)=D(f)(x2),(f(x2)-f(x1))/(x2-x1)=D(f)(x1)}, {x1=0..1,x2=0..1});
x2,x1:=eval([x2,x1], %)[];

plot([f(x),solve((y-f(x1))/(f(x2)-f(x1))=(x-x1)/(x2-x1),y)], x=-0.1..1, color=[red,blue]);

I do not understand the meaning "I want to find the graph of f above in the interval [x1,x2] or [x2,x1] alongside the same function f but which also satisfies at the point of intersection with the above graph f''(x)=0."

Please specify which polynomials are in question. You wrote "How can I define a set generated by two degree n  polynomials with binary coefficients?". But in your example for  n = 2  we see a product of two polynomials, each of the 1st degree and with two variables. And what will happen for n = 3?

Present here the full code of your worksheet.

I think that the point of the problem is that one of the lines is vertical, that is, it is not expressed in the form y = ax + b. If we reflect everything relatively straight line  y=x  (I just changed the order of coordinates), then it is obvious that the angle will not change, but the result will be correct:

restart:
with(geometry):
point(o, 0, 0): 
point(A, 1, 0): 
point(d, 2, 0):
point(F,  1.4472135960, 0.8944271920): 
line(lOD, [o, d]): 
line(lAF, [A, F]): 
alpha := FindAngle(lOD, lAF);

@acer  Thank you for your ways to solving the problem. Using the latter option, I could easily find  R=cos(Pi/17)  in the radicals. Only had to add as an option the field extension  sqrt(17)  in the  factor  command:

restart;
P:=op(1,convert(cos(Pi/17),RootOf));
P1:=factor(P,sqrt(17));
Rs:=simplify(radnormal~([solve(P1,explicit)]),size):
R:=select(r->is(r-cos(Pi/17)=0),Rs)[1];
simplify(R-cos(Pi/17));

@Anthrazit  If you do not want to use custom procedures, then first run  interface(displayprecision=n):  command, where  n  is the number of digits that you want after the decimal point.

Example of use:
interface(displayprecision=3):
evalf(Pi);
                                       3.142

See help on the interface command for details.

@anthonyfl  Very small changes to the code will be required. The animation parameter will be  x  not  y :

restart;
f:=x^(1/2):
g:=x^2/8:
X:=r*cos(phi): Y:=r*sin(phi):
P:=plot3d(eval([[X,Y,f],[X,Y,g]],x=r), r=0..4, phi=0..2*Pi, style=surface, color=["Khaki","LightBlue"], scaling=constrained, axes=normal, labels=[z,x,y], orientation=[20,80], transparency=0.3):
F:=x->plots:-display(plot3d([[X,Y,x^2/8],[X,Y,sqrt(x)]], r=x..x+h, phi=0..2*Pi, style=surface, color=gold), plot3d([eval([X,Y,H],r=x),eval([X,Y,H],r=x+h)], H=x^2/8..sqrt(x), phi=0..2*Pi, style=surface, color=gold)):
h:=0.15:
plots:-animate(F,[x], x=0..4-h, frames=60, background=P);