@Markiyan Hirnyk   Everything works, I checked. On the instrument panel of PP click insert  (not paste)  and so on (see the order of actions above). After  All Files  click your  gif-file from documents.

@laporte bernard   I do not see any errors. To better see the elliptic paraboloid, you just need to extend the domain of definition for the graph of the function  f (x, y)  (the red line on the surface shows the steepest increase of the function  f(x,y) ):

restart;

f:=(x,y)->2+x^2 +1/4*(y^2);

Px, Py:=1, 2;

DPx, DPy:=D[1](f)(Px,Py), D[2](f)(Px,Py);  # Partial derivatives

maximize(DPx*cos(alpha)+DPy*sin(alpha), location);  # Maximum of directional derivative

A:=plot3d(f(x,y), x=-3..3, y=-sqrt(36-4*x^2)..sqrt(36-4*x^2), axes=normal):

B:=plots[arrow]([Px,Py,0],[DPx,DPy,0], color=red):

C1:=plots[spacecurve]( [1,2,t], t=0..4, color=black,linestyle=2):

C2:=plots[spacecurve]( [2*t+1,t+2,f(2*t+1,t+2)], t=0..0.8, color=red, thickness=3):

T:=plots[textplot3d]([[1,2.2,0,"P"], [1,2.2,4,"Q"]], align=right, font=[TIMES,BOLD,16], color=blue):

plots[display](A,B,C1,C2,T, orientation=[-15,65], scaling=constrained);

@Markiyan Hirnyk   Thanks for the clarification of the causes of issues.

@acer  Thanks for the detailed answer! 

@Kitonum  If you replace the fourth line of code with the following line, the procedure will run faster:

L:=[seq(ithprime(k) $ min(n, floor(a/ithprime(k))), k=1..pi(a))];

@fairuzalwani 

E := [0, 2, 7, 15, 26, 40]:

map(C, [3, 10, 24], [0, 2, 7, 15, 26, 40]);

                           [3, 4, 5]

@Carl Love   

E:= [0,2,7,15,26,40]:

1+ListTools:-BinaryPlace(E, Pi);

   Error, (in ListTools:-BinaryPlace) the 2nd argument should be of type {numeric, string}

@one_man   a := remove(i -> i=x+2, a);

@Markiyan Hirnyk   For earlier versions of Maple yours is not appropriate. For example in M 12

seq(RandomTools[Generate](float(range=-1..1, digits=5, method=uniform)), n=1..20);

       Error, (in RandomTools:-Generate) cannot specify a range crossing the origin

So I'm trying to write a version-independent code.

@mahmood180  Yes, I already wrote to you about it.

@mahmood180   Yes, you are right! In the last line of the code replace  M  by  N  (after the sign $) .

What are  а, b  and  c  in your system?  What are the relations  a , b , c  with numbers  a1, a2, a3, a4, a5, a6, a7, a8, a9 ?

@Preben Alsholm   Thank you for your response. But still remains unclear what significance has the singularity for this plotting? Why Maple constructs in this way? After all, for any C>0  the curve  1 / (x ^ 2 + y ^ 2) = C  is a circle.

@royadrb8989   What is  G  and which answer you expect?

@Carl Love   You are right. Easy correction of my code gives the same result:

restart;

n:=0:

S:=combinat[permute]([$ 3..14], 5):

for s in S do

F, H, J, K, L:=op(s):

A := 33 - K - L:

B := 1 + F - J:

C := -15 - F + J + K + L:

d := 15 + H - K:

E := 16 - F - H + J + K:

G := 34 - H - J - L:

i := 18 - J - K:

if convert(map(x->is(x>=3 and x<=14), [A,B,C,d,E,G,i]), `and`)  then n:=n+1: M[n]:=[A,B,C,d,E,G,i,F,H,J,K,L]:  fi:

od:

M:=convert(M,list):

remove(x->nops(convert(x,set))<12, M):

nops(%);  # Total number of solutions

                          0