Multiply function-valued matrices

Question

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[1][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[1][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]&?

Answer 1

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:

Inner[mul,{{a,b},{c,d}},{{w,x},{y,z}},add]

However, Q1 and Q2 are written as functions, rather than matrices with function-valued elements. Using them as written, the matrix product would be implemented like this:

Q = Q1[#1,#2].Q2[#1,#2]&

Then a particular component could be reached with Q[#1,#2][[row,col]]&.