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

Look at these links:

https://www.mapleprimes.com/posts/200188-Congratulations-On-Valentines-Day

https://www.mapleprimes.com/posts/130622-Another-Animation-Of-Heart-Shape

@mahmood1800  Sorry for my inattention. I corrected the definitions of  b[n,m]  in the body of  HybrFunc  procedure. Code1_new  file is replaced by  Code1_new1  (see my answer).

@kivan 

g:=(x,k)->Int(exp(k*cos(z)), z = 0 .. x):
f:=(y,k)->fsolve(y=g(x,k), x=0..infinity):  # The inverse to g  
h:=(u,v,k1,k2)->evalf(Int(exp(sin('f'(x,k1))*sin('f'(y,k2))), [x = 0..u, y = 0..v]));
 

The examples of use:

evalf[15](g(2,3));
f(%,3);
h(4,5,6,7);
                                                       15.2117840729652
                                                           2.000000000
                                                           20.00022605