@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
 

@Annonymouse   Use  seq  command to write this shorter:

S1 := convert~(S, list): 
with(ArrayTools): 
S2 := rhs~(Array(S1));
m, n := op(S2)[1 .. 2]; 
Res:=seq(AllNonZero(S2[i]), i = m);

     
We see that  true  are on the 9th, 10th and 13th places. These 3 numbers we can return programmatically as follows (continuation of the previous code):

ListTools:-SearchAll(true, [Res]);
                                            9, 10, 13

Solve_Problems_MWE2_new.mw

Edit.

@kivan  The reason is that your integral is not expressed in terms of known functions, so it can only be calculated numerically. Therefore, the inverse function can not be expressed symbolically, but can only be calculated numerically using  fsolve  command:

restart;
g:=x->Int(exp(cos(z)), z = 0 .. x): 
h:=y->fsolve(y=g(x), x=0..infinity):   

Example of use:

plot([g,h], 0..10, 0..10, color=[red,blue], scaling=constrained);  # Visualization
evalf[15](g(7));
h(%);

Edit.