To plot, use  plots:-implicitplot  command. Since under the roots there must be positive numbers, your graph will be between 2 hyperbolas  kx*f=sqrt(56.791296)  and  kx*f=sqrt(469.31696) (they are green on the plot). If you want to see a graph for large values of one of the variables, the range of the other variable should be close to  0 .

Download Dispersion_new.mw

You have a problem with boundary conditions. Such problems can very rarely be solved symbolically. For a numerical solution, all parameters must be specified. Also one more condition should be added (I took  theta(0) = 1). The solution is conveniently written in the form of a procedure whose formal parameters are  omega, alpha, M, B, L :

restart;
Sol:=proc(omega,alpha,M,B,L)
local A, BCs, Eq;
A := -(alpha*exp(-sqrt((omega+1)*omega*(M^2+alpha+1))*eta/(omega+1))*omega+alpha*exp(-sqrt((omega+1)*omega*(M^2+alpha+1))*eta/(omega+1))+exp(-sqrt((omega+1)*omega*(M^2+alpha+1))*eta/(omega+1))*omega-alpha*omega+exp(-sqrt((omega+1)*omega*(M^2+alpha+1))*eta/(omega+1))-alpha-omega-1)/sqrt((omega+1)*omega*(M^2+alpha+1));
Eq:=(1+B)*(diff(theta(eta), eta, eta))+C*A*(diff(theta(eta), eta)) = 0;
BCs := (D(theta))(0) = -1,  theta(L)=0, theta(0) = 1;
dsolve({Eq, BCs}, numeric);
end proc:

Example of use:

sol:=Sol(1,2,3,4,5);
plots:-odeplot(sol, [eta,theta(eta)], eta=0..5, scaling=constrained, size=[1000,250]);
 

seq(`if`(sin(k*Pi/4)<>0, sin(k*Pi/4), NULL), k=1..8);

# Or in 2 steps

seq(sin(k*Pi/4), k=1..8);
op(select(`<>`, [%], 0));
 

L:=[1,2,5,4,5,2,5]:
M:=convert(combinat:-choose(L,2), set);
remove(t->t[1]=t[2], M);
nops(%);

# Or

convert(L, set):
n:=nops(%);
n*(n-1)/2;

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]});