@zdamasceno  If I understood correctly, then just add 1 point  [0, (1/4)*Pi-1/2*(b_u-y_f)]:

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:=[(1/4)*Pi-1/2*(b_u-y_f), seq(rhs(fsolve(eval(eqns, f_z=a), s)[]), a=0.00005..0.0005, 0.000001)];
plot([0,seq(a, a=0.00005..0.0005, 0.000001)], S);

@Muhammad Usman  The following code produces the desired result:

restart;

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

for n from 1 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]/sqrt(2)/sqrt(Pi));
od;  od;  od;
A[1];
A[2];
<A[1], Matrix(3); Matrix(3), A[2]>;

Edit.

@waseem  Probably you have not read about what I wrote above: "Unfortunately, the developers did not provide an image of the points in the legends. This can be done manually, but rather cumbersome."

@waseem

Probably you have an older version of Maple, which does not have this option. Try this code:

restart;
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=line, 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]):
          
B:=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=point, color=[red$3,blue$3,black$3],       symbol=[solidbox, solidcircle, soliddiamond], numpoints=9, adaptive=false):

C:=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, C, size=[500,500]);

@vv  Formally, you are right. Because the velocity of motion in this model is equal to the gradient, we obtain a solution that is formally equal to  (0,0)  only for  t=infinity. But for  t=10  we get  (x,y) = (10*exp(-40),10*exp(-20)) = (4.248354255*10^(-17), 2.061153622*10^(-8)) , so that you are already at the top.

It is easy to write equations where the speed is constant, but this complicates the model.

@Muhammad Usman

I did not understand the meaning of this
                                     

@digerdiga  Do you want this output?

ex:=exp(-x^2/2*epsilon)*%sqrt((sqrt(x)/2/h)^2);
                                            

@MarcoFerro 

restart;
evalindets(x^2*y^3+alpha*x*y^5, `*`, p->lcoeff(p, [x,y])*m[degree(p,x),degree(p,y)]);

                                                      alpha*m[1, 5]+m[2, 3]

@MarcoFerro 

Example:

restart;
evalindets(x^2*y^3+x*y^5, `*`, p->m[degree(p,x),degree(p,y)]);

@das1404  Of course I was wrong. There is no second error in your code. However, in principle this statement  "Strings can not act as names" is true, for example:

"x" := 2;
    Error, invalid left hand side of assignment
 

@adxters  This procedure allows us to sort the sublists of some list by increasing the first element in the sublists. See an example:

L:=[[2,1], [1,3], [0,2]]:
sort(L, (a,b)->a[1]<=b[1]);
                                                    [[0, 2], [1, 3], [2, 1]]

@abhilashun  You can easily calculate the volume of this body numerically, using the Monte Carlo method (of course with low accuracy):

f1:=x^4+2*x^2*y^2+2*x^2*z^2+y^4+2*y^2*z^2+z^4-2*x^3-2*x*y^2-2*x*z^2-(79/25)*x^2-(104/25)*y^2-(8/5)*z^2+(104/25)*x:
f2:=x^4+2*x^2*y^2+2*x^2*z^2+y^4+2*y^2*z^2+z^4-2*x^2*y-2*y^3-2*y*z^2-(104/25)*x^2-(79/25)*y^2-(8/5)*z^2+(104/25)*y:
f3:=x^4+2*x^2*y^2+2*x^2*z^2+y^4+2*y^2*z^2+z^4+2*x^3+2*x*y^2+2*x*z^2-(79/25)*x^2-(104/25)*y^2-(8/5)*z^2-(104/25)*x:
f4:=x^4+2*x^2*y^2+2*x^2*z^2+y^4+2*y^2*z^2+z^4+2*x^2*y+2*y^3+2*y*z^2-(104/25)*x^2-(79/25)*y^2-(8/5)*z^2-(104/25)*y:
evalf(Int(`if`(f1<=0 and f2<=0 and f3<=0 and f4<=0, 1, 0), [x=-2..2, y=-2..2, z=-2..2], method = _MonteCarlo, epsilon=0.001));
                                                             12.16037145
 

@abhilashun

plots:-implicitplot3d(max(f1,f2,f3,f4)=0,  x=-2..2, y=-2..2,z=-2..2, style=surface, axes=normal, scaling=constrained, lightmodel=light3, numpoints=1000000); 

@Deank1905  You should more clearly formulate your problem and also download (or give a link) the worksheet.

@acer  I think your approach will be too difficult for a beginner. Because we already know that  f(1)=0 , then we can just write

f := x -> 2*x^2-2;
1/D(f)(1);
                                                

or if we do not know the value of  x  in  y=f(x)  for  y=0 :

f := x -> 2*x^2-2;
x:=eval(x, solve({f(x)=0, x>=0})) ;
1/D(f)(x);