So it sounds like B is a 3xn matrix, where each column {xj,yj,zj} is the components of some separable function fj[x,y,z] = xj + yj + zj?
I'm going to slightly change the definition instead and have a length n vector of functions of x,y,z
{f1[x,y,z], ... , fn[x,y,z]}
I first create a matrix A where a 1 denotes "differentiate with this variable":
vars = {x, y, z};
nVars = Length@vars;
A = With[{id = IdentityMatrix[nVars]},
id~Join~(id + RotateLeft[id])];
A // MatrixForm
Then calculate the Jacobian and transpose. Here for example with 2 functions:
n = 2;
funs = Table[Symbol["f" <> ToString[i]] @@ vars, {i, n}];
J = Outer[D, funs, vars] // Transpose;
J//MatrixForm
Then dot A with the transposed Jacobian J:
(A . J) // MatrixForm
We can alternatively do this without naming the variables, and just specifying which argument we want to differentiate with using A:
ClearAll["`*"]
nVars = 3;
id = IdentityMatrix[nVars];
A = id~Join~(id + RotateLeft[id]);
n = 2;
funNames = Symbol["f" <> ToString[#]] & /@ Range[n];
J = Outer[Derivative[Sequence @@ #1][#2] &, id, funNames, 1, 1];
A . J // MatrixForm
There is likely a much cleaner way to do this using Inner, but I haven't figured out a nice way of doing it yet.
