Thank you for the answer.
So, the reason for opposite orders will remain a mistery.

Happy New Year! (actually also for all the mapleprimes members).

V.A.

@Markiyan Hirnyk 

I don't think that sort could be useful for an answer.
I'd like to know whether there is a consistent/logical rendering order for the objects in display.
The CURVES seem to be on top of POLYGONS and the orders differ.
When using transparency e.g., the order is important.
 

@torabi 

I suspect that the system has no solution (even after trying some other initial conditions).
Unfortunately it is hard to prove that this is the case and even harder to find an existence theorem for such a system. Sorry.
(If the system came from a concrete problem, you should try anoher approach).

@Rouben Rostamian  

Thank you Rouben, I'll keep these worksheets in my "special" collection.
In order to solve the memory problem for animation, one may use Explore. This way a single frame is present in memory but now the animation speed depends on  the computer speed.

anim := proc(tau, lambda, T)
  local R := 6;  # orbit's radius
      display([
      sphere([0,0,0], 1, style=surface, color="Orange"),
      tubeplot([R*cos(t), R*sin(t), 0], t=0..2*Pi, radius=0.04, color=red, style=surface),
      translate(globe(tau, lambda, T), R*cos(T), R*sin(T), 0)],
      view=[-R-1..R+1,-R-1..R+1,-1..1],
      scaling=constrained, lightmodel=light4, axes=none, orientation=[-144,65,0]
    )
end proc:
Explore(anim(tau,lambda,T),
parameters=[[tau=0..2*Pi,animate=false,shown=true], [lambda=0..2*Pi,animate=false,shown=true], [T=0..2*Pi,animate=true]],
initialvalues=[tau=Pi/6,lambda=Pi/4], loop, size=[1000,1000], numframes=240, autorun);

@torabi 

So, what is your conclusion?

I seem to remember that you have posted this system several times, without success.
Are you sure that it is not similar to a simpler one:

dsolve({diff(f(x),x)=1/g(x), diff(g(x),x)=0,  g(0)=0, f(0)=0});

?

@Rouben Rostamian  

Very nice solution, vote up!

PS. An animated 3d plot would be also nice; maybe someone will have the patience to do it.

@one man 
I thought that you want to know how your surfaces really look (globally).
But if you want to test the procedure then please choose examples which satisfy the mentioned limitations.

@one man 
For these implicit surfaces the above procedure cannot be used (for a global representation) because the conditions are not satisfied (about interior point and star-shapeness).
But it is not very difficult to parametrize them.
1.   f:=(x^2+y^2-2/5)^2+(z+sin(x*y+z))^4-1/10;

plot3d([(1/10)*(40+10*10^(1/2)*cos(u))^(1/2)*cos(v),
(1/10)*(40+10*10^(1/2)*cos(u))^(1/2)*sin(v), 
-(1/100)*(40+10*10^(1/2)*cos(u))*cos(v)*sin(v)+
RootOf(-cos(u)*10^(1/2)*cos(v)*sin(v)-10^(1/2)*(10^(1/2)*abs(sin(u)))^(1/2)-4*cos(v)*sin(v)+10*sin(_Z)+10*_Z)],
u=0..2*Pi, v=0..2*Pi, orientation=[60,40]);

2. f := x+y+z(+const);
Probably it's a joke.

3. f:= z-1/2 * exp (sin (x + 5/2 * y + z));

Y:=(k,x,z) -> -(2/5)*x-(2/5)*z+(2/5)*( (-1)^k*arcsin(ln(2*z)) + k*Pi);
plot3d( [seq(Y(k,x,z),k=0..5)], x=0..4*Pi, z=exp(-1)/2 .. exp(1)/2, labels=["x","z","y"],
         scaling=constrained, orientation=[-75,4,-114]);

# The graph has an infinity of connected components, each one being unbounded

@Ian Thompson 
But here c is not a table; c[0], c[1],... are simply indexed names.

@Markiyan Hirnyk 
This is obviously a bug. The smart algorithm in implicitplot is supposed to detect curves and not to add points not satisfying the equation. In this case it inserts in the plot some other curves which have nothing to do with the equation.

@Axel Vogt 

Download a0b0c0.mw

@Axel Vogt 

If f is homogeneous (actually only positively), any solution is a positive multiple of a solution from [-1,1]^3. I hope that this is clear.

Now, indeed, for c=-1 there exist solutions with a,b outside [-1,1], but we are not interested in these. They will be automatically included in the general solution. You can check easily that for some k and u (in the answer above) your solution will be obtained [you did not provide the values for a,b, otherwise I could show you the values for k and u and the permutation].

@Carl Love 
- The reason for which I did not try to use evalhf or other optimizations is that I think that anyway fsolve is the bottleneck.
- The equation is not linear; it is  f(X0 + r*V) = 0 where X0, V are constants in R^3.
- I am not sure if the body is star-shaped. It would not be complicated to test this (numerically) inside the procedure: fsolve should be replaced with a procedure which gives a warning if there are more than one root in 0..R, but probably this would be too time-consumming.

Edit: I tested with an adhoc fsolve and the body seems to be star-shaped.

@Ronan 
If you look at the exact solution (see my answer below), it should be clear that Maple (or any other CAS) has no chance to find it by itself.