L:=[1,2,5,4,5]:
m:=max(L);
map(p->`if`(p=m, 0, p), L);

Your example:

solve(sin(x)=y, x, allsolutions);
                                            -2*arcsin(y)*_B1+Pi*_B1+2*Pi*_Z1+arcsin(y)

The parameters:  _Z1  is an integer,  _B1  is equal  0  or  1  .

Use  fnormal  command for this (and also  abs  command to delete minus sign before zero if necessary):

abs(fnormal(eval(f, x=ans[1])));
abs(fnormal(eval(f, x=ans[2])));
                                                         0.
                                                         0.

See help on  fnormal  command for details.
 

We can verify this identity for specific values of parameter  s  using  identify command.

Example:

restart;
s:=5:
evalf(product(ithprime(j)^s/(ithprime(j)^s-1), j=1..1000));
identify(%);
                                           

Addition. Sure this is not proof from the standpoint of pure mathematics, but from the point of view of common sense it is quite convincing.

Let  A  is your matrix. Try

simplify~(A) assuming positive;

factor(expand(numer(z)))/denom(z);
                                                                 x^a

# Or at once

factor(expand(z));
                                                                 x^a
 

In order for Maple to solve the system numerically, the older derivatives must be explicitly expressed through the lower derivatives and the remaining variables. Maple can easily cope with this if the system is linear with respect to the higher derivatives. In the third equation we have a quadratic equation with respect to the highest derivative. We can solve this equation and take one of the solutions, but then a new problem arises - for your parameter values, in the process of solving negative numbers appear in the power of 1/2. See the attached file.

New_new.mw

To verify this, we first find  Re(f(x)) , then the denominator of it, and then check in two ways that on the interval  0.3..0.36  it takes only positive values, that is, it is nowhere equal to zero.

Addition. On the plot, you can clearly see this minimum:

plot(A, x=0.344..0.345, 5*10^244..7*10^244, numpoints=10000);

Download Asymptotes_new.mw

restart;
p:=dsolve({D(y)(x)=(-2*exp(x))/(2+exp(x)), y(0)=1/9}, numeric);
plots:-odeplot(p, [x,y(x)], x=0..20);

This equation can be solved symbolically (exactly):
q:=dsolve({D(y)(x)=(-2*exp(x))/(2+exp(x)), y(0)=1/9});

Use  eval  command  or  D  differential operator for this. The latter is preferable:

diff(x^2*y^2+3*x*y^2, x, y);
eval(%, [x=1,y=2]);
f:=(x,y)->x^2*y^2+3*x*y^2;
D[1,2](f)(1,2);
 

I guess that matrices  A[1]  and  A[2]  should be generated of the numbers  Omega[m, j]  . Then we get

restart; 
alpha := 1: k := 2: M := 3:
printlevel := 3;

for n while n <= 2^(k-1) do

for m from 0 while m <= M-1 do

for j from 0 while j <= M-1 do

Omega[m, j] := 2^((1/2)*k)*sqrt(GAMMA(j+1)*(j+alpha)*GAMMA(alpha)^2/(Pi*2^(1-2*alpha)*GAMMA(j+2*alpha)))*(sum((-1)^i*GAMMA(j-i+alpha)*2^(j-2*i)*(sum((1/2)*binomial(m, l)*(2*n-1)^(m-l)*(1+(-1)^(j-2*i+l))*GAMMA((1/2)*j-i+(1/2)*l+1/2)*GAMMA(alpha+1/2)/GAMMA(alpha+1+(1/2)*j-i+(1/2)*l), l = 0 .. m))/(GAMMA(alpha)*factorial(i)*factorial(j-2*i)), i = 0 .. floor((1/2)*j)))/2^(k*(m+1));
A[n]:=Matrix(3, (i,j)->Omega[i-1,j-1]);
od;  od;  od;
A[1];
A[2];
Matrix(2, 2, {(1, 1) = A[1], (1, 2) = 0, (2, 1) = 0, (2, 2) = A[2]});

max[index]  command appeared in one of the recent versions of Maple. If you have an older version, you can use other methods. Below for your example, 2 more ways:

restart;
L:=[1,4,7,4,2,1]: m:=max(L):

# The first way
for i from 1 while L[i]<m do  
od:
L[1..i];

# The second way
L[1..ListTools:-Search(m, L)];

Use  plot  and  discont commands for this.

Example:
plot(sinh(x), x=-3..3);
discont(sinh(x), x);  # The set of points of discontinuity

s:=t->piecewise(t>=0, sin(Pi*t), 0):
g:=t->piecewise(s(t)>=1/2, 1, 0):
plot([s(t), g(t)], t=0..3.5, color=[red,blue], thickness=3, symbol=solidcircle, discont=true, scaling=constrained, size=[1000,400], legend=["s(t","g(t)"]);

You can add  style=pointline  option to your code:

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, style=pointline, linestyle = [solid,longdash,dashdot],'
          thickness = 2',color=[red$3,blue$3,black$3],
          legend=["Nb=0.1",""$2,"Nb=0.2",""$2,"Nb=0.3",""$2], symbol=[solidbox, solidcircle, soliddiamond], numpoints=9, adaptive=false):
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, size=[500,500]);

Addition. Unfortunately, the developers did not provide an image of the points in the legends. This can be done manually, but rather cumbersome.