Multiply function-valued matrices

I would like to create a notebook for calculations in supersymmetric quantum mechanics. The basic building blocks are two functions (or maybe better: operators):

A := (W[#1] #2 + D[#2, #1]) &;
B := (W[#1] #2 -D[#2, #1]) &;

With these definitions you can apply these operators to a pair of a variable x and a function f[x] like A[x,f[x]] and one obtains f[x] W[x] + Derivative[f][x].

Now, these functions A and B are embedded into matrices:

Q1 := ({{0, 0}, {A[#1, #2], 0}}) &;
Q2 := ({{0, B[#1, #2]}, {0, 0}}) &;

Again, it is possible to apply these matrices to a variable x and a function f[x] like Q[x,f[x]] and one obtains

{{0, 0}, {f[x] W[x] + Derivative[f][x], 0}}

or in traditional form I wonder, if there is a possibility to multiply the matrices Q1 and Q2 in such a way that the result is again a function that expects the two arguments [x,f[x]].

The trouble is, that the internal matrix operations have to distinguish between multiplication and addition of ordinary numbers and functions. I succeeded to implement addition and multiplication that is sensitive to the nature of the operands on the level of A and B,eg.

mul[a_, b_] := If[NumberQ[a] && NumberQ[b], a*b,
If[NumberQ[a], Through[a*b[#1, #2], 0] &,
If[NumberQ[b], Through[b*a[#1, #2], 0] &,
a[#1, b[#1, #2]] &
]]];

and

add[a_, b_] :=If[NumberQ[a] && NumberQ[b], a + b,
If[NumberQ[a], a + b[#1, #2] &,
If[NumberQ[b], b + a[#1, #2] &,
Through[a[#1, #2] + b[#1, #2], 0] &]]]
]];

I have no idea, how to deal with the matrix multiplication. For instance, how can I get the (2,1) component of Q as a function. If I write Q[[1,2,1]], I get #2 W[#1] but I would need W[#1] #2 + D[#2,#1]&?

• A = W[#1] #2 + HoldForm@D[#2, #1] & ? – Fortsaint Jun 17 at 21:20

If you're looking for a generalized matrix product, then Inner will do the trick. This snippet performs a matrix product, where multiplication is replaced with mul and addition with add: