It's a brilliant Answer, Markiyan. I vote up. Here are two slight improvements:

1. The Asker wanted 60 contours.
2. By including a 1% expansion multiplier in the transform, we insure that the contour lines remain visible as solid lines above the surface.

with(plots): with(plottools):
f:= (psi-Pi/2)^2 + (theta-.5)^2;   
a:= plot3d(1, psi= 0..2*Pi, theta= 0..Pi, coords= spherical,
     grid= [77$2], style= patchnogrid, color= f
);
b:= contourplot(f, psi = 0..2*Pi, theta= 0..Pi, contours= 60, thickness= 0);
g:= transform(
     (psi, theta)-> 
          1.01*[(cos,sin)(psi)*~sin(theta), cos(theta)]
):
display([a, g(b)], axes= frame, scaling= constrained, lightmodel= NONE); 

It's a brilliant Answer, Markiyan. I vote up. Here are two slight improvements:

1. The Asker wanted 60 contours.
2. By including a 1% expansion multiplier in the transform, we insure that the contour lines remain visible as solid lines above the surface.

with(plots): with(plottools):
f:= (psi-Pi/2)^2 + (theta-.5)^2;   
a:= plot3d(1, psi= 0..2*Pi, theta= 0..Pi, coords= spherical,
     grid= [77$2], style= patchnogrid, color= f
);
b:= contourplot(f, psi = 0..2*Pi, theta= 0..Pi, contours= 60, thickness= 0);
g:= transform(
     (psi, theta)-> 
          1.01*[(cos,sin)(psi)*~sin(theta), cos(theta)]
):
display([a, g(b)], axes= frame, scaling= constrained, lightmodel= NONE); 

@Markiyan Hirnyk 

This must be a bug in DEtools[autonomous] such that it does not recognize the k and m as constants. Consider

restart:
constants:= constants, k, m:
sys:= {
     diff(x(t),t) = -x(t) + 2*z(t)*y(t) - 2*k*z(t)^2,
     diff(y(t),t) = -y(t) + k*m - z(t)*(x(t)/2 - 1),
     diff(z(t),t) = m
}:
DEtools[autonomous](sys, {x,y,z}, t);
                              true

@Markiyan Hirnyk 

This must be a bug in DEtools[autonomous] such that it does not recognize the k and m as constants. Consider

restart:
constants:= constants, k, m:
sys:= {
     diff(x(t),t) = -x(t) + 2*z(t)*y(t) - 2*k*z(t)^2,
     diff(y(t),t) = -y(t) + k*m - z(t)*(x(t)/2 - 1),
     diff(z(t),t) = m
}:
DEtools[autonomous](sys, {x,y,z}, t);
                              true

I think that you have the roles of theta and phi switched. Try changing the plot command to

plot3d(f, theta= 0..2*Pi, phi= 0..Pi, coords= spherical);

Are these the plots that you expected?

There is no attached file.

In your worksheet, you use both N and n. Did you intend for those to be the same? Does U have columns 0..m? or 1..m? or something else?

@pallav Sorry that I did not anticipate this already. Try this:

restart:
a := 2; b := .29; d := 1.85;
k:= 0;
for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     k:= k+1;
     tau[k] := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
end do:
convert(tau, list)[];

@pallav Sorry that I did not anticipate this already. Try this:

restart:
a := 2; b := .29; d := 1.85;
k:= 0;
for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     k:= k+1;
     tau[k] := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
end do:
convert(tau, list)[];

@pallav Then you need to save the taus while in the loop, and use a statement after the end of the loop to display the values. Like this:

restart:
a := 2; b := .29; d := 1.85;
for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     tau[h] := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
end do:
convert(tau, list)[];

@pallav Then you need to save the taus while in the loop, and use a statement after the end of the loop to display the values. Like this:

restart:
a := 2; b := .29; d := 1.85;
for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     tau[h] := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
end do:
convert(tau, list)[];

Could you give an example? It is not clear to me why you presented the lists three times.

If they are lists, you should enclose them in []: [a,a,a,b,c]. Otherwise, they are called sequences, which are more difficult to work with than true lists.

Do you want the rules to apply to sublists also? That is, should [b,a,a,a,b,c,e] become [b,b,b,b,c,e,e]?

@Markiyan Hirnyk That is not an example, and you know it! I ask for a solution set returned by solve which contains a trivial equation. And that is what the Asker is referring to in his followup to your Answer.

@Markiyan Hirnyk That is not an example, and you know it! I ask for a solution set returned by solve which contains a trivial equation. And that is what the Asker is referring to in his followup to your Answer.

@Markiyan Hirnyk Please kindly show me an example of a solution set returned by solve, or a similar command, which contains a trivial equation and yet represents only a finite number of solutions.