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

@Markiyan Hirnyk 
OK, I'll watch this thread waiting to see a solution on your taste.

The parametrization I have obtained looks a little better but I am still not satisfied and I think that it would be very difficult to find an acceptable general method. The discretization (mesh) is too "odd". The picture is in the worksheet, you can rotate it an see the irregular grids.

-param3d.mw

@one man 

AFAIK the method is about natural parametrization of implicit curves in R^(n+1), a curve being given as an intersection of n implicit surfaces. You have applied it for n=2.
Do you have a reference for the original paper?
PS. Why do you spell Draghilev? The author is probably Russian. But the spelling AFAIK is not unique in such cases, just think about Chebyshev/Tchebysheff/Cebâşev. And in Russian, all the western names are also spelled phonetically!
[In Zentralblatt MATH he appears as "Dragilev, A.V."]

@one man 
It is not that complicated. The metod is the same as the one I have posted above except that you obtain first a discrete parametrization of the implicit surface (you have a post about this).

P.S. I wrote 2 procedures inspired by your method. Actually the main ingredient is the Dragilev method for the natural parametrization of an implicit curve in R^3.
It works OK for local parametrizations (obtaining a MESH structure) but for global surfaces it is hard to obtain a convenient initial curve and even harder  to find the "orthogonals" other than the planar ones you have used. (The planar "orthogonals" are not good enough; you may convince yourself trying to parametrize a rotated ellipsoid.)
That's why I have not posted them.

BTW, to color both sides of an implicit surface one may use e.g.

f1:=(x1^2+x2^2-0.4)^2+(x3+sin(x1*x2+x3))^4-0.1:
a:=plots:-implicitplot3d([f1=0], x1 = -1 .. 1., x2 = -1 .. 1., x3 = -0.1 .. 1.9, color=[green], numpoints=10000, style=surface):
b:=plots:-implicitplot3d([f1=0.001], x1 = -1 .. 1., x2 = -1 .. 1., x3 = -0.1 .. 1.9, color=[red], numpoints=10000, style=surface):
plots:-display(a,b,orientation=[37,61,-171]);