In tensor calculation, I need to do the following thing:

Define a derivative operator `Drv[fun,i]`, such that

 Drv[f[i],j]=f[i,j]
 Drv[f[i,j],k]=f[i,j,k]
 Drv[f[i]**g[j],k]=f[i,k]**g[j]+f[i]**g[j,k]


----------
**UPDATE**

I think this question in fact is quite important to me, and the answer, however is not quite complete, thus I will try to make the question more clearly.

The basic setting is that we want to define a *derivative* rule for the operator `NonCommutativeMultiply`, recall what's a derivative, call `Drv`:

 1. To do the derivation, we must know who is/isn't a function with respect to the derivative variable, since `Drv` acts on them differently. Thus as first step, it should be declare a set of functions that are real functions, let's say `f, g, h`.
 2. The basic rule for a derivative is that:

 - Linearity: `Drv[c**f]=c**Drv[f]` and `Drv[f+g]=Drv[f]+Drv[g]`;
 - Distributive: `Drv[f**g]=Drv[f]**g+f**Drv[g]`;

 3. for multivariable functions, we just write $Drv[f(x_1,x_2,...,x_n),x_k]$ as $f[k]$, and $Drv[Drv[f(x_1,x_2,...,x_n),x_i],x_j]$ as $f[i,j]$ and so on. For example:
$$
 Drv[f,i]=f[i]\\
 Drv[f[i],j]=f[i,j]\\
 Drv[f[i,j],k]=f[i,j,k]
$$
the properties of derivative is just reads:

**linearity**
$$
 Drv[c**f,i]=c**Drv[f,i]=c**f[i], 
 Drv[f[i]+g,j]=f[i,j]+g[j]
$$
for c is a function which is independent on $x_i$. But when $c$ is a real number, 
$$
 Drv[2 f,i]=2f[i]
$$
**distributive**
$$
Drv[f**g[i]**h[j,k],l]=f[l]**g[i]**h[j,k]+f**g[i,l]**h[j,k]+f**g[i]**h[j,k,l].
$$