Curves:=plot([tan(x), tan(Pi/4)+D(tan)(Pi/4)*(x-Pi/4), arctan(x), arctan(1)+D(arctan)(1)*(x-1)], x=0..3, 0..3, color=[red,yellow,blue,green], discont, legend=[tan(x),"Tangent for tan(x)", arctan(x),"Tangent for arctan(x)"]):
Points:=plot([[[Pi/4,1]], [[1,Pi/4]]], style=point, color=[red,blue], symbol=solidcircle, symbolsize=10):
plots:-display(Curves, Points, scaling=constrained, size=[400,400]);

BVP := [4*(diff(u(x, t), t))-9*(diff(u(x, t), x, x))-5*u(x, t) = 0, u(0, t) = 0, u(6, t) = 0, u(x, 0) = sin((1/6)*Pi*x)^2]:
sol:=pdsolve(BVP):        
U := unapply(rhs(sol), [x, t]);  
plot3d(add(op(1,U(x, t)), n=1..20, 2), x = 0 .. 6, t = 0 .. 4);

To find the derivative of the function  f  with respect to the first variable for certain values of these variables, you must first define the function  f  itself in terms of its variables, for example:

restart;
f:=(u,v)->u^2+u*v;
D[1](f)(z*sqrt(a/nu[f]), U*t/(2*x));

# The first way
seq(seq(exp(a*Pi/b), b=1..10), a=1..15);
nops([%]); 

# The second way
k:=0:
for a from 1 to 15 do
for b from 1 to 10 do
k:=k+1; L[k]:=exp(a*Pi/b);
od: od:
convert(L, list)[];  

The results are identical.

Place:=proc(L,n)
local k;
uses ListTools;
k:=Search(n,L);
if k<>0 then return k else NULL fi;
end proc:

Examples of use:

Place([1,2,4,5], 4);
Place([1,2,4,5], 6);

 

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;