Product between operators - Unwanted behavior when replacing operator by matrices

[Michele Cotrufo]

I define two operators and their product

from sympy.physics.quantum.operator import Operator

A,B = Operator('A'), Operator('B')

C=A*B

later, I want to replace the operators A and B by two squared matrices of the same dimensions. In particular, I want C to be equal to the product of the matrices associated to A and B.
However, it seems that sympy does something different:

Asub = sp.Matrix([[0,1],[1,2]])

Bsub = sp.Matrix([[1,1],[1,3]])

C.subs([(A,Asub)])

this last line returns 

Matrix([

[0,   B],

[B, 2*B]])

so it considered B as a scalar, and it multiplied each element of the matrix A by B

If instead I try

C.subs([(A,Asub),(B,Bsub) ])

I get  (sorry for the bad formatting)

Matrix([

[                        0, Matrix([

[1, 1],

[1, 3]])],

[Matrix([

[1, 1],

[1, 3]]), Matrix([

[2, 2],

[2, 6]])]])

so essentially it returns a 2x2 BLOCK matricx.
How can I force sympy to treat A*B as a matrix product?