I suggest using the functional style instead  eval : the code turns out to be more compact. This is also convenient if you want to calculate the values of  lambda  function at individual points. I also added text labels for your groups of graphs as in your picture:

h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
K1:=(4/h(z)^4)-(sin(alpha)/F)-h(z)^2+Nb*h(z)^4:
lambda:=(F,Nb,delta2)->Int(K1,z=0..1):

L1:=0.2:
d1:=0.2:
alpha:=Pi/6:
A:=plot( [seq(seq(lambda(F,Nb,delta2), Nb=[0.1,0.2,0.3]), F=[0.1,0.2,0.5])], delta2=0.02..0.1, color=[red$3,blue$3,black$3]):
B:=plots:-textplot([[0.05,-1,F=0.1],[0.05,1.5,F=0.2], [0.05,3,F=0.5]], font=[times, 14]):

plots:-display(A, B);

The output:

plot([4*t, 4*t-0.5*t^2, t=0..10]);

See  wiki  and the help on  ?plot,details

findOrderOf:=proc(g,p)
  local i:=1, gpwr:=1;
  for i from 1 to p-1 do
    gpwr:=gpwr*(g mod p);
    if gpwr = 1 then return i
    end if;
  end do;
end proc:

primes:=select(isprime,[`$`(2000 .. 3000)]);
select(p->findOrderOf(2,p)=p-1 and findOrderOf(3,p)=p-1,primes);

This works:

Download test_(1)_new.mw

In the body of the procedure  poly_out  replace the line

plots[display](trigr, lbase, trapezr, outAr, axes=none, scaling=constrained);

by the line

plots[display](trigr, lbase, trapezr, outAr, axes=normal, scaling=constrained);

We see that Maple reduces the calculation of this integral to a definite integral, which probably can not be calculated symbolically, but it is easily can be done numerically. Here is this calculation, as well as visualization: the cylindrical body is plotted, the volume of which is the computed integral:

Download Double_Integral.mw

restart;
q:=Matrix(3,10,(i,j)->i^2+j+2*i); 
with(plots):
V:=i->(q[..,round(i)]):
animate(arrow,['V'(k), width=[0.4,relative=false],head_width=[1,relative=false],color=red], k = 0.5..10.45, frames=60, axes=normal, paraminfo=false);
 

solve({7*cos(2*t)=7*combine(cos(t)^2)-5, t>=0, t<2*Pi}, t, allsolutions, explicit);  # Symbolic
evalf(%);  # Numeric
                 

Let the function  f  be the inverse of the function temperature vs chemicals, M  is your Array. Do the following:

A:=convert(M, Matrix);
L:=convert(<f~(A[..,2]) | A[..,1]>, listlist);
plot(L);
 

eq16:=r(t)=d[vol]*V/(KUS*V^2+L*tau);
r(t)=(numer(rhs(eq16))/V)/'(expand(denom(rhs(eq16))/V))':
eq17:=%;

Slightly reducing the number of digits, it is easy to check the identity of the lists:

evalb(evalf[8](listA) = evalf[8](listB));
 

restart;
f:=unapply(<r(t)*cos(t), r(t)*sin(t)>, t):
v:=diff~(f(t), t);  # The vector of velocity

Example of use. The calculation of the vector of velocity for the curve  r(t)=t  (Archimedes' spiral) at the point  t=1 :

f:=unapply(<t*cos(t),t*sin(t)>, t):
v:=unapply(diff~(f(t),t), t):
v(t);  # The vector of velocity
v(1);  # Velocity at the point  t=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), scaling=constrained);  # Visualization

Edit.

Obviously, this is a bug. Here's a very trivial example:

f1 := x->x^2:   f2 := x->x^2:
is(f1=f2);
                                                false

Probably  is  command simply does not work with procedures. Obviously to check the equality of the two procedures  f  and  g , it suffices to verify that  f(x)=g(x)  for an arbitrary  x:

f1:=x->x^2:  f2:=t->t^2:
is(f1(a) = f2(a));
                                                 true

Use  evalc  command for this.

Example:

evalc(exp(k*Pi*I));
                                       cos(k*Pi)+I*sin(k*Pi)

For some reason, D  operator does not work in  Physics  package, but  diff  command is working:

restart:
with(Physics[Vectors]):
Setup(mathematicalnotation = true):
r1:=t->2*t^2*_i+16*_j+(10*t-12)*_k;
r2:=t->(20-6*t)*_i+4*t^2*_j+2*t^2*_k;
v1:= unapply(diff(r1(t),t), t);
v2:= unapply(diff(r2(t),t), t);
arccos(v1(2).v2(2)/(Norm(v1(2))*Norm(v1(2))));  # Angle in radians
evalf(%*180/Pi);  # Angle in degrees