@asa12 
Let us try to formulate mathematically the problem (if I understand it correctly).

You have a function space  FS, for example C^infinity(J)  (= the algebra of infinitely differentiable functions on the open interval J) . You also have an operator A defined (partially) on FS, in your case it is:

   A(a) = a' / (a * a' - a).

You want to find another operator f defined on FS such that

   f( A(a) + A(b) )  = f(A(a)) * f(A(b)),  for each a,b in FS.

Such an f is for example f = exp, or more generally, f(a) = exp(s*a),  s being a scalar.

Now you want an operator f, similar to diff. I suppose that this means:

  f( A(a) * A(b) ) = f(A(a))*A(b) + A(a)*f(A(b)).      [*]

(I suppose you want to exclude the trivial case f=0).

Maple probably cannot do this. And I am not sure if such an operator - other than f = D - exists. It may depend on FS, which you should know, because it is your ploblem. And I am not sure if [*] is what you want.

@asa12 

Maybe you want an Infinitesimal generator, see https://en.wikipedia.org/wiki/C0-semigroup#Infinitesimal_generator
But you must learn no ask the question in proper terms (preferably using some maths and without non-working Maple code), otherwise it is "quasi-impossible" to be understood.
 

@Carl Love 

Nice solution, vote up.
Don't you thing that using multi-indices is more efficient?
I wonder why the OP has no reaction.

@Kitonum 

It's OK but you should add the step about the equivalence, for those thinking that it is not obvious.

[ 2y=x(1-y^2)  <==> 2y/(1-y^2) = x,  because y cannot be +-1]

P.S. In my opinion, using Maple without knowing maths is not recommended.

@Markiyan Hirnyk 

Of course you do, it's not something new.

@Markiyan Hirnyk 

They include the cases when both sides are undefined (but actually solutions in the Riemann sphere).
If you want to exclude this cases, simply take integer _Z
with  - _Z <> W + 2*W^2,  for integer W.

@one man 

For a cone the geodesic distance is straightforward, so the exact parametric equation of the spiral is obvious.

@Carl Love 

Looking at showstat(`mod/Factors`)  we can see that a table s is created and then
factors := convert(s,'multiset');
As shown in my answer (and reply), this will (could) produce a random order.

To test this, I'd suggest to modify `convert/multiset` inserting indexorder and see whether the random behavior persists.

By coloring the cube with color f(x,y,z) you will have information on f only on the boundary of the cube.
That is, only f(0,y,z), f(1,y,z),...,f(x,y,1) are visualized.
You can't "see" e,g. the value f(0.5, 0.6, 0.7).
So, you will see the max/min only if they are on the boundary.

@Carl Love 

Sorry for the confusion, but I thought that the reference to Eigenvectors is obvious because of the presence of lambda.

@Carl Love 

Eigenvectors produces a random order for the eigenvalues because after solving the characteristic equation with the local variable lambda, indices is called.
You are right about indexorder. The procedure which calls indices without indexorder is `convert/multiset`. If we change it inserting indexorder, the random behavior disappears.

@Markiyan Hirnyk 

OK, it's enough for me too.
P.S. I wish you good luck in obtaining small Groebner bases using the cited articles.

@Markiyan Hirnyk 

I was not saying that the paper is not serious. But their approach is only for (some) numeric problems. The Groebner bases are used for much more!

@Markiyan Hirnyk 

I did not read them because as I have added (probably you did not see because it was during your reply) I am interested in general Groebner bases which are inherently huge.

Edit. OK, you made me curious and I have looked over [1]. As I have anticipated, they approximate a Groebner basis using Newton's method. That is not what most people want!

@Markiyan Hirnyk 

OK, but if for 70 equations we need supercomputers then ...

Edit. I was (and am) interested in general Groebner bases, not in some versions designed for numerics. Here the situation is unfortunately clear: the Grobner basis is unique and very often is HUGE.