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

But the angle at x=a is not  arctan(|f'(a) - g'(a)|). It should be  | arctan(f'(a)) - arctan(g'(a)) |.

PS. I doubt that a beginner who writes cosx^2  will grasp map2, ~, @

@Carl Love 

@wyoum 

Try:

subs( Physics:-`*` = :-`.`, expr );

subs( Physics:-`.` = :-`.`, expr );

@wyoum 
If you are interested in FP groups, the Physics package will be anyway not enough.
With GroupTheory you will be able to simplify much better.
Example:

with(GroupTheory):
G := < a, b | a^2 = 1, b^4 = 1, (a.b)^2 = 1 >:
E:=Subgroup([[]],G);  # trivial subgroup
Factor( a.b.a.b^3 .a.b.a, E);
             [1,b]

So, in the group G one has:    a.b.a.b^3 .a.b.a = b

I think I have made an error while writing my answer. Trying to detete it and reload the worksheet, I think I have deleted  Markiyan Hirnyk' s answer.
Sorry, Markiyan, please repost the answer.