It may be the order of evaluation:

coor = {x, y, z};
list = {#1, #2, #3};
e[i_, j_] :> ((D[Evaluate[u[i] @@ list], coor[[j]]] + 
      D[Evaluate[u[j] @@ list], coor[[i]]])/2 &)

you can see that coor and list are not inserted into the rule. If you use With the values will be injected:

With[{coor = {x, y, z},
  list = {#1, #2, #3}},
 e[i_, j_] :> ((D[Evaluate[u[i] @@ list], coor[[j]]] + 
       D[Evaluate[u[j] @@ list], coor[[i]]])/2 &)]

Another approach:

sost[n_Integer] := 
  e[i_, j_] :> 
   Function[
    Evaluate[((Derivative @@ UnitVector[n, j])[u[i]] @@ 
         Array[Slot, n] + (Derivative @@ UnitVector[n, i])[u[j]] @@ 
         Array[Slot, n])/2]];
e[1, 2][x, y, z] /. sost[3]
D[e[1, 2][x, y, z], y] /. sost[3]

Dear Glanluca, First thanks again for the reply. Sorry I didn't make it clear in my topic. But the whole point of setting up two lists "coor={x,y,z} and list={#1,#2,#3}" is to generalize to arbitrary dimensions, as described in the original post. The 3D test expression I give is just a special case.

I am trying to figure out why my substitution doesn't work - when e(i,j) is substituted by the derivative of a specified function u[i][x,y,z], it worked. When e(i,j) is substituted by (Evaluate[u[i] @@ list] &), it worked. But when I try to substitute e(i,j) with derivative of Evaluate[u[i] @@ list], it doesn't work. Why is that? Is there any way to do it correctly?

Thank you!

Hi Gianluca, Thank you for your reply. I have tried your solution on this test sample:

sost = e[i, j] :> Function[((Derivative @@ UnitVector[3, j])[ u[i]][##] + (Derivative @@ UnitVector[3, i])[u[j]][##])/2]; D[e[1, 2][x, y, z],y] /. sost

Which returns 0 on my computer. So I suppose it doesn't work.

For the second solution, I agree it will work out if I enumerate all possible e[i,j] like that, but is there any way to avoid enumeration using the substitution rules?

Thanks!