@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. 

@toandhsp  For convenience, I collected all the codes in one file.

Triangles.mw

@waseem  In the code below, I took into account your wishes:

restart:
h1:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
h2:=z->1-(delta2/2)*(1 + cos(2*(Pi/L2)*(z - d2 - L2))):
h3:=z->1+(delta2/2):
assign(seq(K[i]=((4/h||i(z)^4)-(sin(alpha)/F)-h||i(z)^2+Nb*h||i(z)^4), i=1..3)):

lambda1:=Int(K[1],z=0..0.2):
lambda2:=Int(K[2],z=0.2..0.4):
lambda3:=Int(K[3],z=0.4..0.6):
lambda:=lambda1+lambda2+lambda3:
 
F:=0.3:
L1:=0.2:
d1:=0.2:
d2:=0.2:
L2:=0.3:
alpha:=Pi/6:
plot( [seq(eval(lambda, Nb=j), j in [0.1,0.2,0.3])], delta2=0.02..0.1);

@Markiyan Hirnyk  It is obvious that both results (mine and vv's) are correct, but relate to slightly different distributions.

@waseem 

You made several syntax errors:

1. It is not useful to assign one name  h  to different objects, because only the last assignment will work.

2. This line of code causes an error

h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))): for 0<z<0.2

Just take  for 0<z<0.2  away. This range is indicated in the integral.

@Markiyan Hirnyk  In my answer I did not use this function, so I did not look at Maple for help on it and immediately looked at the wiki.

@Markiyan Hirnyk  Quotation from the wiki:

"In probability theory and statistics, the geometric distribution is either of two discrete probability distributions:

• The probability distribution of the number X of Bernoulli trials needed to get one success, supported on the set { 1, 2, 3, ...}
• The probability distribution of the number Y = X − 1 of failures before the first success, supported on the set { 0, 1, 2, 3, ... }"      You did not specify what kind of case is meant. I took the first as the more common .

@escorpsy I think that it is impossible to obtain analytically explicit dependences of roots on parameters, but it is possible to investigate how the roots change when individual parameters change. This can be done, for example, using  Explore  command. See update to my answer.