@zdamasceno  Now everything is clear. At first I just did not understand the meaning of what you first called the boundary condition. In fact, this is just an initial approximation using Newton's method. In Maple, you do not need to specify it. Here are 2 solutions to your example. In the first, this initial approximation is indicated. In the second solution, the range for the wanted solution is specified instead. The results (lists  S1  and  S2) are identical.

Download eqns.mw

@zdamasceno  For any value of the parameter  f_z  , you get a trigonometric equation for  s - variable, which has the infinite number of solutions. The condition   s= (1/4)*Pi-1/2*(b_u-y_f)  when f_z=0  allows us to choose some branch of solutions. But if we substitute  f_z=0  into the equation, we get an error (the divizion by 0). I guess that you have incorrectly set  Z.  f_z - variable must be in the numerator (not in the denominator). In this case we have  s= (1/4)*Pi-1/2*(b_u-y_f)=0.3783306793   for  f_z=0 :

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

@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);