@waseem   Change the last line to

plots:-display(Vector([seq(plot([seq(eval(lambda(F), Nb=j), j=[0.1,0.2,0.3])], delta2=0.02..0.1, legend=[Nb=0.1,Nb=0.2,Nb=0.3], title=('F'=F)), F=[10,20,30])]), size=[900,500]);

@ganelon  From this specification nothing will change: the solution is unique for any  y>0

@acer  How did you know about [Gr,Br] ?

@Carl Love  You're right.

@das1404   If I need to draw several graphs on one plot, then I specify the common options (axes, scaling, view) only at the end of  plots[display]  command. In this case there will be no contradiction between the options.

@Preben Alsholm 

We can simplify a little more by

EXPR2:=combine(EXPR1, ln) assuming h>0, R0>0;
length(EXPR2);

@alfarunner  There are no problems here. The coefficients before rho hat and phi hat are the coordinates of the vector v(t) in this basis (here phi=t):

restart;
f:=unapply(<r(t)*cos(t),r(t)*sin(t)>, t):
v:=unapply(diff~(f(t),t), t):
V:=<cos(t),sin(t); -sin(t),cos(t)>.v(t);
simplify(V);  # The coordinates in the new basis (rho hat and phi hat)

I have never worked in the Physics package.

@alfarunner  Above I wrote how to calculate the coordinates of the velocity vector in standard Cartesian coordinates. If you need the same, but in other coordinates, you just need to do a coordinate transformation, multiplying by the corresponding matrix:

restart;
f:=unapply(<r(t)*cos(t),r(t)*sin(t)>, t):
v:=unapply(diff~(f(t),t), t):
V:=<cos(t),sin(t); -sin(t),cos(t)>.v(t);

f:=unapply(<t*cos(t),t*sin(t)>, t):  # The same example
v:=unapply(diff~(f(t),t), t):
v(t);
v(1);
simplify(<cos(1),sin(1); -sin(1),cos(1)>.v(1));
plots:-display(plot([convert(f(t),list)[],t=0..3], color=red, thickness=3), plots:-arrow(f(1),v(1), width=[0.01,relative=false],head_width=[0.08,relative=false],color=blue), plots:-arrow(f(1),<cos(1),sin(1)>, width=[0.01,relative=false],head_width=[0.08,relative=false],color=black), plots:-arrow(f(1),<-sin(1),cos(1)>, width=[0.01,relative=false],head_width=[0.08,relative=false],color=black), scaling=constrained);
 

@Preben Alsholm  Thank you for clarifying. Yet from the standpoint of pure mathematics it seems that this is not a very good design.

@minhthien2016  To do this, the code should be rewritten as a procedure. After that, you can apply this procedure to any cube by specifying as a procedure argument the list of cube vertices.

Here is this procedure:

restart;
RegTetr:=proc(L)
local T, Dist, P;
uses combinat;
T:=combinat:-choose(L, 4);
Dist:=(p,q)->sqrt((p[1]-q[1])^2+(p[2]-q[2])^2+(p[3]-q[3])^2);
P:=v->combinat:-choose(v, 2);
[seq(`if`(nops(convert((Dist@op)~(P(t)),set))=1, t, NULL), t=T)];  
end proc:

Example of use:

L:=[[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1], [1,0,1], [1,1,1], [0,1,1]]:
L1:=2*~(map(t->t-~[1/2,1/2,1/2], L));
RegTetr(L1);

@vv  Your solution is impressive, but weakly motivated (much remains behind the scenes). Why did you take these ranges  x=50..100, y=30..100, H=400..500 ? In addition, we can not so round up  y=40.00000000, because we obtain 0 under the logarithm sign.

Here is another solution in which we first eliminate  H , simplify the obtained system , and then decide, taking as the initial point x=41, y=41, because obviously should be  x>40, y>40  for real solution :

restart;
Digits:=1000:
Sys:=[21000 = ((10/1000)/60)*4181*983*(x-y),
H = -((10/1000)/60)*4181*ln(1-((x-y)/(x-40))),
H = 124+102*ln(x-40)+(7+2*ln(x-40))];
A:=eliminate(Sys, H):
Sys1:=convert(A[2], list):
Sys2:=expand(Sys1) assuming x>40, y>40;
fsolve(Sys2, {x=41,y=41}):
Sol:={%[], eval(Sys[3],%)};  # The final result
eval(Sys, Sol):  # Check
evalf[20](%);
 

Download System.mw

@Aaeru Michi  Maple 16 just does not cope with this example (I solved in Maple 2017). Here are 2 ways that maybe help you:

minimize(V, x = 0 .. 1, location) assuming x>0;  # or

minimize(-(1/2)*abs(x)*(c*x^2-2*L+x)/(2*c*x+1)^2, x = 0 .. 1, location);  # I just removed  signum(2*c*x+1)  from your example

@Aaeru Michi   From help: "The Optimization package is a collection of commands for numerically solving optimization problems, which involve finding the minimum or maximum of an objective function possibly subject to constraints.  The package takes advantage of built-in library routines provided by the Numerical Algorithms Group (NAG)"

See help on  ?Optimization  for details.

@fereydoon_shekofte 

Look at this:

A := plottools[polygon]([[sqrt(3)+1, 0, 0], [1, 1, 0], [1, -1, 0]], color = "Violet"):

B := plottools[polygon]([[sqrt(3)+1, 0, 0.01], [1, 1, 0.01], [1, -1, 0.01]], color = "Green"):

These are two different polygons at the distance = 0.01  one from the other.

@exality  I never use the context menu and I always work by direct commands in 1D input mode. It's more convenient, faster and more reliable.