@emma hassan I hope that you caught the small mistake in my Answer, which I've since corrected. I had

S:= LinearSolve(A, XY);

which should have been

Solution:= LinearSolve(A, XY);

What is mdcs? I guess it's a debugger.

What is mdcs? I guess it's a debugger.

I got somewhat different results. In Maple 17 Standard (Windows8, 64bit), I had no trouble: The help pages open and look great. In Maple 16.02 Standard (Windows8, 64bit), I get "No matches found" for both ?Janet and help(Janet), even though the package commands themselves are loaded.

I don't have Matlab to test this, but what if you take out the quotes? Make it

maple(fsolve(a*x^2-4))

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)[];