You can use  plots:-dualaxisplot  command for this:

restart;
X:=[10000,70000,75000,80000,85000,90000,95000]:
Pi_b:=[12478250.8,18963457.2,17846277.6,16743622.7,15655492.6,14581887.3,13522806.7]:
Pi_r:=[484495.57,204495.57,240745.572,281162.238,325745.572,374495.572,427412.238]:
plots:-dualaxisplot(plot(X,Pi_b, style=pointline,  color=red, labels = [X, pi__b^`1`], symbol=solidcircle), plot(X,Pi_r, style=pointline, color = blue, labels = [X, pi__r^`1`], symbol=solidcircle), size=[900,500]);

Consider a cube with a side of 7 and divide its 3 faces into rectangles 1 x 7, 6 x 7, 4 x 7, 3 x 7, 5 x 7, 2 x 7. These rectangles are the bases of 6 pyramids. Connect the vertices of each rectangle with a single vertex of the cube. We get 6 pyramids and it is easy to check that they satisfy the conditions of the problem. In fact, each of these pyramids has one side of the base equal to 7 and the height equal to 7. Therefore, their volumes are related as third sides of the bases, that is, 1 : 2 : 3 : 4 : 5 : 6   :

times := [-2.019619977*Unit('s'), 2.019619977*Unit('s')]:
select(t->coeff(t, Unit('s'))>0, times)[];

When you find intersection points, for example, of  MF1  line with the ellipse, you get  2  points, one of which will be point  M . But it will not necessarily be the first of the two solutions returned. Therefore, from the two solutions  sol , you must choose the point that does not coincide with point  M .

restart;
Fig := proc(t) 
local a, b, c, courbe, sol, sol1, dr, tx;
uses geometry, plots; 
_EnvHorizontalName := 'x'; _EnvVerticalName := 'y'; 
a := 11; b := 7; c := sqrt(a^2 - b^2); 
ellipse(e1, x^2/a^2 + y^2/b^2 = 1); 
point(Oo, 0, 0); point(M, a*cos(t), b*sin(t)); 
point(F1, -c, 0); point(F2, c, 0); 
line(MF1, [M, F1]); 
line(MF2, [M, F2]);
 
sol:=solve({Equation(MF1), x^2/a^2 + y^2/b^2 = 1}, {x, y}, explicit);
point(H1,eval([x,y],sol[1])[]); 
if evalf(distance(H1,M))<10^(-6) then point(H1,eval([x,y],sol[2])[]) fi;  
line(MH1, [M, H1]); 

sol:=solve({Equation(MF2), x^2/a^2 + y^2/b^2 = 1}, {x, y}, explicit); 
point(H2,eval([x,y],sol[1])[]); 
if evalf(distance(H2,M))<10^(-6) then point(H2,eval([x,y],sol[2])[]) fi;  
line(MH2, [M, H2]); 

courbe := implicitplot(x^2/a^2 + (a^2 + c^2)^2*y^2/b^2 - 1 = 0, x = -a .. a, y = -b .. b, color = cyan); 
tx := textplot([[coordinates(M)[], "M"], [coordinates(Oo)[], "O"], [coordinates(H1)[], "H1"], [coordinates(H2)[], "H2"], [coordinates(F1)[], "F1"], [coordinates(F2)[], "F2"]], font = [times, bold, 16], align = [above, left]); 
dr := draw([e1(color = blue), MH1(color = magenta), MH2(color = magenta), M(color = red, symbol = solidcircle, symbolsize = 12), H1(color = red, symbol = solidcircle, symbolsize = 12), H2(color = red, symbol = solidcircle, symbolsize = 12), F1(color = red, symbol = solidcircle, symbolsize = 12), F2(color = red, symbol = solidcircle, symbolsize = 12), Oo(color = red, symbol = solidcircle, symbolsize = 12)]); plots:-display([dr, tx, courbe], scaling = constrained, axes = normal, title = "Ellipse et normales ", titlefont = [HELVETICA, 14]); 
end proc:

plots:-animate(Fig, [t], t = 0.1 .. 2*Pi, frames = 150);

Edit.

The updated last line of your code:

P := [eval([k1, delta][], s[2]), s[1]]; 
display(plot3d(TRC(k1, delta), k1 = 0 .. 1, delta = 0 .. 1), pointplot3d(P, color = red, symbol = solidcircle, symbolsize = 17), textplot3d([P[], "P"], font = [times, 16], align = {ABOVE, LEFT}), orientation = [45, 75], axes = normal);

 Do  labels=["δ0", s__2]

You wrote in the worksheet : "As shown in Section 2 all of the centres of the STCs lie on a circle, Dd+, related to Dd. "  This statement is incorrect. I took a large circle of radius 4 and the next circle of radius 3 with center C1. Direct calculation found 3 more circles with centers C2, C3, C4. Then found: a circle circumscribed around points C1, C2, C3 (red color) and a circle circumscribed around points C2, C3, C4 (blue color). We see that these circles are different.

restart;
local O:
R:=4.: r1:=3.:
dist:=(A,B)->sqrt((A[1]-B[1])^2+(A[2]-B[2])^2):
center:=(A,B,C)->solve({(B-A)[1]*(x-(A+B)[1]/2)+(B-A)[2]*(y-(A+B)[2]/2),(C-A)[1]*(x-(A+C)[1]/2)+(C-A)[2]*(y-(A+C)[2]/2)}, {x,y}):
O:=[0,R]: O1:=[x1,y1]: O2:=[x2,y2]: O3:=[x3,y3]: 
O1:=eval(O1,solve({dist(O,O1)=R+r1,y1=r1}, {x1,y1},explicit)[1]);
O2:=eval([x2,y2],simplify(solve({dist(O1,O2)=r1+r2, dist(O,O2)=R+r2,y2=r2},{x2,r2,y2}, explicit)[1])):
O2:=[2*sqrt(R)*sqrt(r1)*sqrt(R)/(sqrt(R)+sqrt(r1)),r1*(sqrt(R)/(sqrt(R)+sqrt(r1)))^2];
k:=simplify(sqrt(R)/(sqrt(R)+sqrt(O2[2]))) assuming R>0,r1>0:
O3:=simplify([O2[1]*k,O2[2]*k^2]) assuming R>0,r1>0;
k:=simplify(sqrt(R)/(sqrt(R)+sqrt(O3[2]))) assuming R>0,r1>0:
O4:=[O3[1]*k,O3[2]*k^2];
c1:=simplify(eval([x,y],center(O1,O2,O3)));
c2:=simplify(eval([x,y],center(O2,O3,O4)));
with(plottools): with(plots):
assign(seq(P||i=circle(O||i,(O||i)[2]), i=1..4)):
C1:=circle(c1,dist(c1,O1), color=red):
C2:=circle(c2,dist(c2,O2), color=blue):
display(circle(O,R),P1,P2,P3,P4,C1,C2, scaling=constrained, view=[0..12,0..10], size=[800,800]);

circles1.mw

PS.  It seems that the centers of all smaller circles lie on the parabola  y=x^2/16  (for the given values  R=4, r1=3):

display(circle(O,R),P1,P2,P3,P4,plot(x^2/16, x=0..10, color=green), pointplot([O1,O2,O3,O4],color=red,symbol=solidcircle,symbolsize=7),scaling=constrained, view=[0..10,0..9], size=[800,800]);

The  evalc  command performs calculations on complex-valued expressions assuming that all parameters are real-valued (here  u(x,t)  and  v(x,t) ) :

restart;
alias(u = u(x, t), v = v(x, t)):
q:=u+I*v:
Eq:=evalc(diff(q, t)+diff((diff(q, x))/evalc((1+abs(q)^2)^(3/2)), x$2)=0);

evalc(Re(lhs(Eq)));
evalc(Im(lhs(Eq)));

Sum  is an inert sum, and to get the numeric result you can use the  value  command:

restart;
n1 := 122:
x := 54:
n2 := 65:
y := 14:
alpha1 := 1.1:
beta1 := 1.1:
alpha2 := 1.1:
beta2 := 1.1:

Sigma2 := (j, l, psi0) -> Sum(Sum(Sum(psi0^(-y - alpha2 - l - alpha2*h - v)*pochhammer(beta2 - 1, h)*pochhammer(beta1 - 1, u)*pochhammer(y + alpha2 + l + alpha2*h, v)*2^(-y - alpha2 - l - alpha2*h - x - alpha1 - j - alpha1*u - v)*(psi0 - 1)^v/(h!*(y + alpha2 + l + alpha2*h)*u!*v!*(y + alpha2 + l + alpha2*h + x + alpha1 + j + alpha1*u + v)), u = 0 .. upto), h = 0 .. upto), v = 0 .. upto):

F2 := psi0 -> Sum(binomial(n1 - x, j)*Sum(binomial(n2 - y, l)*(-1)^(j + l)*Sigma2(j, l, psi0), l = 0 .. n2 - y), j = 0 .. n1 - x):

upto:=10:

value(F2(5.620));

This can be done in many ways. In such cases, I usually write the simplest procedures using standard formulas:

restart;
dist:=(A,B)->sqrt((A[1]-B[1])^2+(A[2]-B[2])^2):
midp:=(A,B)->(A+B)/2:

# Example
dmf1:=[-3, 1]:  dmf2:=[3, 2]:
dist(dmf1,dmf2);
midp(dmf1,dmf2);

If we introduce coordinates  C(0,0), D(486,0), A(x1,y1), B(x2,y2) , we get an example of an underdetermined problem, since we have 4 unknowns, but only 3 conditions (3 known distances) to find them. Such systems usually have an infinite number of solutions. To find a specific solution, it is necessary to assign a certain value to one unknown. I considered 2 cases: x2=120  and  x2=140 . To find the area, I used the Area procedure, a detailed review of it is available at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-

restart;
local D:

Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
add(P[i], i = 1 .. nops(L)); 
end proc:

C:=[0,0]: D:=[486,0]: A:=[x1,y1]: B:=[x2,y2]:
ds:=(M,N)->(M[1]-N[1])^2+(M[2]-N[2])^2:
eqs:=ds(C,B)=(162.8)^2, ds(B,A)=(71.8)^2, ds(A,D)=325^2;
sol1:=fsolve({eqs,x2=120}, {x1=0..infinity,y1=60..infinity,y2=0..infinity,x2=0..infinity});
sol2:=fsolve({eqs,x2=140}, {x1=0..infinity,y1=60..infinity,y2=60..infinity,x2=0..infinity});
L:=[[C,D,A,B,C]]:

Area(eval(L, sol1)); # First solution
Area(eval(L, sol2));  # Second solution

See https://math.stackexchange.com/questions/513071/tetrahedron-inequality?newreg=b9e402d644c3454d92b79b0dd7d32a25

Below in the code, for brevity, I limited myself to the first 10 prime numbers:

restart;
with(combinat):
primes := [seq(ithprime(i), i = 1 .. 10)];
S:=choose(primes, 6):
F:=(x,y,z,u,v,w)->`if`(LinearAlgebra:-Determinant(<0,1,1,1,1; 1,0,x^2,y^2,z^2; 1,x^2,0,u^2,v^2; 1,y^2,u^2,0,w^2; 1,z^2,v^2,w^2,0>)>0,true,false): k:=0:
for s in S do
for t in [seq(seq([s1[],s2[]], s2=permute(convert(s,set) minus convert(s1,set))), s1=choose(s, 3))] do
x,y,z,u,v,w:=t[];
if F(t[]) and min(x+y+z,x+v+u,w+y+u,w+v+z)>max(x+w,y+v,z+u) then k:=k+1; L[k]:=t; break fi;
od: od:
L:=convert(L, list); 
nops(L);

We received the list  L  of 22 solutions.

I've seen many cases where this command fails to work if the system contains inequalities. The workaround is obvious:

restart;
L:=[isolve(a*b=4)];
select(s->`and`(op(rhs~(s)>=~1)), L)[];

PS. For such cases, the  IntegerPoints2   procedure from the post  https://mapleprimes.com/posts/202542-Integer-Points-In-Curved-Regions is useful:

IntegerPoints2({a*b=4,a>=1,b>=1},[a,b]);

@nm  I see that you are more interested in plotting integral curves (stream lines) than just arrows. Below I show how to plot 10 curves for your example. It is not very difficult to plot arrows along each line as well, or just arrows, using the  plots:-arrow  command. The picture clearly shows that the origin is a stable focus.

restart;
with(plots):
sys:=[diff(x(t),t) = -3*x(t)-4*y(t), diff(y(t),t) = 2*x(t)+y(t)]:
S1:=seq(dsolve({sys[],x(0)=4,y(0)=i}), i=0..4):
S2:=seq(dsolve({sys[],x(0)=-4,y(0)=i}), i=-4..0):
P:=plot([seq(eval([x(t),y(t),t=0..2*Pi],s), s=[S1,S2])], x=-4.5..4.5, y=-5..5, color=red, thickness=2):
V1:=seq(seq(arrow(eval(eval([x(t),y(t)],S1[k]),t=T),eval(eval([diff(x(t),t),diff(y(t),t)],S1[k]),t=T),length=0.9,head_width=0.15,head_length=0.2,color=red,border=false), T=[0.05,0.15,0.3,0.6,1.2]), k=1..5):
V2:=seq(seq(arrow(eval(eval([x(t),y(t)],S2[k]),t=T),eval(eval([diff(x(t),t),diff(y(t),t)],S2[k]),t=T),length=0.9, head_width=0.15,head_length=0.2,color=red,border=false), T=[0.05,0.15,0.3,0.6,1.2]), k=1..5):
display(V1,V2,P, size=[800,800], labels=[``,``], font=[times,16]);

The answer has been significantly edited.

stream_lines.mw

The  legend  option is not supported when drawing arrows. Instead, you can make text labels near each vector. In my opinion, this looks even better than using legends:
 

Download legend_question_new.mw