I wonder what are the steps that led to the set of indices and the set of values. And I wonder what you plan on doing with the resulting set of equations. I suspect that you could be about to attempt something terribly inefficient that could be better handled with Matrices.

@OffshoreEngineer I assume from your dependent variable names that you are using cylindrical coordinates. How about trying to solve the BVP in Cartesian coordinates?

@OffshoreEngineer For what it's worth, I tried two other things, with no substantial change in the results. The first was to reduce each coefficient to a single significant digit. The second was to divide through each equation by its leading coefficient. So, I don't think that the phenomenon under discussion in this thread is a numerical/computational issue.

@OffshoreEngineer For what it's worth, I tried two other things, with no substantial change in the results. The first was to reduce each coefficient to a single significant digit. The second was to divide through each equation by its leading coefficient. So, I don't think that the phenomenon under discussion in this thread is a numerical/computational issue.

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

@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?