You can use  textplot  command instead of  legend  option to place the legends in any desired place. Of course, this will require a few more efforts.

Example:

with(plots):
A:=plot([sin(x), cos(x)], x=-Pi..3*Pi, color=[red,blue]):
B:=textplot([[3.7,0.8,sin(x)], [-1.8,0.5,cos(x)]], font=[times,bold,
14]):
C:=textplot([[3,1,__, color=red], [-2.5,0.7,__,color=blue]], font=[times,bold,24]):
display(A, B, C, size=[950,200], scaling=constrained);

Use  PolyhedralSets  package. See the help on this package for details. Here is an example of the plotting of your first set (all inequalities must be non-strict):

with(PolyhedralSets):
P:=PolyhedralSet({40 <= y, x <= y-40, z <= y-40});
Plot(P, color=khaki, axes=normal, view=[-30..60, 0..90, -40..50], orientation=[40,70]); 

We see that this is an unbounded trihedral body.
 

A lot of syntax errors. Also note that  implicitplot  command only works for two variables. To examine the dependence of the plot on the third variable (I took  p  as the third variable), you can use  animate  or Explore  commands:

restart; 
with(plots):
omega[0]:=1: epsilon[1]:=1.04493: epsilon[2]:=0.93259: delta:=0.02:
eq:=(omega[0]*a-a*Omega^2+(3/4)*epsilon[2]*a^3)^2+(delta*Omega*a)^2 = p^2;
animate(implicitplot, [eq, Omega = 0.5 .. 2, a = 0 .. 5, color = blue, thickness = 2, gridrefine = 3], p=0..0.2, frames=120);

# Or
Explore(implicitplot(eq, Omega = 0.5 .. 2, a = 0 .. 5, color = blue, thickness = 2, gridrefine = 3), parameters=[p=0...1.], numframes =120);

Edit.

We can noticeably improve the perception if we emphasize the ellipse by its border and/or a color.  fieldstrength=log   option is also useful, which makes it possible to reduce the difference in the size of the arrows:

restart; 
with(plots):    
f := (x,y) -> `if`(x^2+y^2/2 <= 1, x^2+y^2, 0):     
g := (x,y) -> `if`(x^2+y^2/2 <= 1, x+y, 0):      
plots:-display(plots:-implicitplot(x^2+y^2/2-1<=0, x=-1.5..1.5,y=-1.5..1.5, thickness=0, filledregions=true, coloring=["LightYellow",white]), fieldplot( [f,g], -1.5..1.5,-1.5..1.5, arrows=SLIM, color="Red", fieldstrength=log));

Your algorithm may not provide a minimum spanning tree. You should at each step look for the nearest vertex not for the last vertex of the already constructed tree, but to search for the nearest vertex by scanning all the vertices of the already constructed tree. This will be Prim's algorithm  https://en.wikipedia.org/wiki/Prim%27s_algorithm

Here is an example where your algorithm fails (start from vertex a):

A:=plot([[0,0],[-1,-1],[0,2],[1,3]], style=point, symbol=solidcircle, color=red, symbolsize=20):
B:=plots:-textplot([[0,0,a],[-1,-1,b],[0,2,c],[1,3,d]], color=blue, align={below,right}, font=[times,bold,20]):
plots:-display(A, B, scaling=constrained, axes=none);

Addition. Prim's algorithm is easy to implement in Maple, but now I do not have time to do this (maybe tonight I'll do it)

To improve accuracy, I used timestep and spacestep options. I took the range  t=0..0.5, because already for  t = 0.5  your function is practically equal to 0 (see below). For integration, a spline of 1 order was used :

restart; 
k:=5:
EQ:=diff(u(x,t),t)=k*diff(u(x,t), x$2):
ibc:=u(0,t)=0, u(1,t)=0, u(x,0) = x:
sol:=pdsolve({EQ}, {ibc}, numeric, timestep = 0.001, spacestep = 0.005):
p1:=sol:-plot(u, x=0.5, t=0...0.5, style = line, color = "Blue", legend = "heat Plot", axes=boxed);
M:=op(1, op(1, p1)):
f:=CurveFitting:-Spline(M, t, degree=1):
int(f, t=0..0.5);  # The area under p1
eval(f, t=0.5);  # The value of the function for t=0.5

extrema  command by  's'-option  allows you to find all the critical points (in which both partial derivatives are 0). Student[MultivariateCalculus][SecondDerivativeTest]  command allows you to classify these points (LocaMin, LocalMax, or Saddle) :

restart;

f:=(x,y)->x^3 + 3*x*y^2 - 15*x + y^3 - 15*y:
extrema(f(x,y), { }, {x,y}, 's');  # The minimum and maximum values of the function at critical points

s;  # All critical points

for p in s do
p1:=convert(p, list);
print(Student[MultivariateCalculus][SecondDerivativeTest](f(x,y), lhs~(p1)=rhs~(p1)));
od:  # Classification of the critical points

with(plots):
display(plot3d(f(x,y), x=-6..6, y=-6..6, view=-60..60, orientation=[145,60]), pointplot3d(map(t->[rhs~(t)[ ], f(rhs~(t)[ ])]
, s), symbol=solidcircle, color=red, symbolsize=20));  # Plot

Edit.

f:=x->Statistics:-PDF(Weibull(6.2, 2), x):
A:=plot(f, 0..26, color="DodgerBlue", thickness=3, tickmarks=[spacing(5), spacing(0.02)], axis[2]=[gridlines=8], view=[0..30,0..0.16], size=[600,400], legend=Series1, legendstyle=[font=[times,roman,16], location=right]):
B:=plot([seq([x,f(x)], x=0..26)], style=point, symbol=soliddiamond, symbolsize=18, color="DodgerBlue"):
plots:-display(A, B);

Suppose that such a function  F  exists. Then we have for  t=exp(1) :

F(F(exp(t))) = F(t) = F(exp(1))=1  instead of  0 .

We see that these conditions (F(exp(t)) = t  and  F(F(exp(t))) = 0)

contradict one another.

y=0.0000125698-0.0000125698*cos(54*x);
a:=convert(%, fraction);
b:=content(rhs(a), cos(54*x));
c:=primpart(rhs(a), cos(54*x));
y=b*``(c);
convert(%, float, 6);
 

You just have to solve equation  g(x)=1.2*f(x) :

f:=x->75000*exp(-0.045*x)+100000:
g:=x->300000/(1+100000/3333*exp(-0.08*x)):
fsolve(g(x)=1.2*f(x), x=0..70);  # The result
plot([f, g], 0..70, color=[red,green], legend=[urban,rural]);  # Visualization
 

Try

int(evalc(exp(I*x*(2-y))),[x=-infinity..infinity, y=0..2]);  # Or
evalf(Int(evalc(exp(I*x*(2-y))),[x=-infinity..infinity, y=0..2]));

The  F  procedure solves the problem:

F:=(i, x)->rand(i-x .. i+x)( );

Example of use:

seq(F(20, 5), k=1..20);

            19, 15, 16, 22, 24, 16, 19, 15, 17, 23, 16, 16, 17, 19, 18, 18, 19, 25, 22, 25

Since there are  i+x-(i-x-1)=2*x+1  numbers in total in the range  i-x .. i+x, the probability of each individual number is  1/(2*x+1) .

Edit.

Use initial conditions and  numeric  option. You can get the values of a function at individual points and plot its graph.

Example:

Sol:=dsolve({-1.699678499+diff(z(t), t, t)-3.813322790*z(t)+.6966019659*z(t)^2+0.1548004368e-1/z(t)^2 = 0, z(0)=1, D(z)(0)=0}, z(t), numeric);
Sol(1);
plots:-odeplot(Sol, [t,z(t)], t=0..10, view=0..10);

Obviously this is bug. Even if you use  RealDomain:-solve  command, the error persists :

RealDomain:-solve(-sqrt(k*x) = sqrt(k*x+k-1), {k,x});

                               {k = 1, x = x}