1. The 1000-th prime is  7919  not  1979

2. Use  unapply  command instead of  -> :

restart;
P:= 7919; # (this is the1000-th prime)
S:=5000; k:=5:
a:=[0$i=1..k];
f:=[0$i=1..k];
a[1]:=S;
f:=unapply(a[1], x);

with(RandomTools):
for i from 2 to k do
a[i]:=Generate(integer(range=1..P-1));
f:=unapply(f(x)+a[i]*x^(k-1), x);
print(f(x));
od;

restart;
sys := {6*(diff(a(t), t))^2+12*a(t)*(diff(a(t), t$2))-3*a(t)^2*phi(t)^(-2*c)*sqrt(1-alpha*(diff(phi(t), t))^2), 2*c*a(t)^3*phi(t)^(-2*c-1)*sqrt(1-alpha*(diff(phi(t), t))^2)-3*alpha*a(t)^2*phi(t)^(-2*c)*(diff(a(t), t))*(diff(phi(t), t))/sqrt(1-alpha*(diff(phi(t), t))^2)-alpha*a(t)^3*phi(t)^(-2*c)*(diff(phi(t), t$2))/sqrt(1-alpha*(diff(phi(t), t))^2)+2*c*alpha*a(t)^3*phi(t)^(-2*c-1)*(diff(phi(t), t))^2/sqrt(1-alpha*(diff(phi(t), t))^2)-alpha^2*a(t)^3*phi(t)^(-2*c)*(diff(phi(t), t))^2*(diff(phi(t), t$2))/(1-alpha*(diff(phi(t), t))^2)^(3/2)};
R(t) := 6*((diff(a(t), t))^2/a(t)^2+(diff(a(t), t$2))/a(t));
W(t) := -phi(t)^(-2*c)*sqrt(1-alpha*(diff(phi(t), t))^2)/(1/a(t)^3+a(t)^3+phi(t)^(-2*c)/sqrt(1-alpha*(diff(phi(t), t))^2));
inc:=a(0)=1.5,D(a)(0)=0.1,phi(0)=1,D(phi)(0)=-0.5;
sol:=dsolve({sys[], inc}, numeric, parameters=[c, alpha]);

sol(parameters=[1,1]);
plots:-odeplot(sol, [[t,a(t)], [t,R(t)]], t=0..10, color=[red,blue]);
 

I guess that  :=  (assignment sign) should be instead of  =  :

restart; 
M := 5;
for i while i <= M do N[i, 0](u) := 1; N[0, i](u) := 0 end do;

Another way is the usage  assign  command for a multiple assignment (without any loops):

 assign(seq(op([N[i, 0](u)=1, N[0, i](u)=0]), i=1..M));
 

restart;
f:=(x,c,k)->c*x/(k+x):
plot([f(x,5,2), f(x,5,3)], x=0..50, color=[red,blue], legend=["c=5 \n k=2","c=5 \n k=3"], size=[800,400]);

                       

The corrected version:

1_(1)_new.mw

restart;
eqns := {(1-sin(b_u)*sin(s)/(cos(b_u-y_f)*cos(s-y_f)))*(1/cos(s-y_f)^2-1/sin(s)^2) = -(cot(s)+tan(s-y_f)+Z)*sin(b_u)*(cos(s)/cos(s-y_f)+sin(s)*sin(s-y_f)/cos(s-y_f)^2)/cos(b_u-y_f)};
b_u := 1/tan(0.8);
Z := 892/(27417000*f_z);
y_f := 9*Pi/180;
S:=[seq(rhs(fsolve(eval(eqns, f_z=a), s)[]), a=0.00005..0.0005, 0.000001)];

S is the list of desired solutions.

Addition. The plot:

plot([seq(a, a=0.00005..0.0005, 0.000001)], S);

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;