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

@Carl Love 
Yes, I was once the author of such an "infringement", but I cannot understant how posting the result of a Maple command (showstat in this case) could be interpreted as a copyright problem.

@Carl Love 

plots:-transform works if the surface is simple wrt one of the axes (y in this case).
It will not work e.g. for
S1:= plot3d([1, theta, phi], theta= 0..3*Pi/2, phi= 0..Pi, coords= spherical):
But in such cases, for general parametric plots a small normal perturbation will work
(see http://www.mapleprimes.com/posts/203796-Equidistant-Surface-)

CP:=(v,w)->[v[2]*w[3]-v[3]*w[2], v[3]*w[1]-v[1]*w[3], v[1]*w[2]-v[2]*w[1]]: #Cross Product
IP:=(v,w)->v[1]*w[1]+v[2]*w[2]+v[3]*w[3]: # Inner Product
UV:=v -> v/~sqrt( IP(v,v)): #Unit Vector

r:=[sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)]:  # Example
Nr:=UV( CP( diff(r,theta), diff(r,phi) ) ):  # Unit Normal vector to r
Er:=simplify(r + 1/200*Nr):;           # Equidistant to r

plot3d( [r,Er], theta= 0..3*Pi/2,  phi= 0.001..Pi,color=[red,green] );

The procedure `hypergeom/check_parameters` checks whether the upper parameters which are integer <=0 are at least as many as the lower ones. [maybe this could be relaxed].

I have noticed that:
limit( hypergeom([1, -1, 1/2], [k,-3], 1),  k=-12 ); 
      71/72

@lyakhovda 
But my previous reply answers to your question for r=23 too; it is NO.

@John Fredsted 

The OP wants to choose himself the free parameters.
Actually for a random linear system, this can be done with probability 1, but not for each and every system.

@lyakhovda 

If the rank is 6, then there will be 29-6=23 free variables, so {Z20,...,Z28} are not enough.
And you cannot know whether e.g. Z24 can be a free variable because the system could contain the equation Z24 = 7.

Yes, Maple has big problems with oscillating integrals.

J := int(abs(cos(1/t)), t = 0 .. 1);
is wrong (undefined).
And it is a challenge to compute J with 10 (or 30) digits in Maple!