Couldn't you simply substitute X and Y into the polynomial P as the values of x and y? (Throwing worries about ordering to the wind because we have that they commute.)
sage: # set things up
sage: R.<x,y> = QQ[]
sage: P = R.random_element(10,12)
sage: M = MatrixSpace(QQ, 5)
sage: X = M.random_element()
sage: Y = M.random_element()
sage:
sage: # check the inputs
sage: P
-x^3*y^7 + 1/2*x^2*y^6 + 1/2*x^7 - 13*x^5 - 1/2*y^5 + 2*x*y^2 + y^3
sage: X
[ 1/2 0 -1 0 -1/2]
[ -2 0 -1/2 1 -1/2]
[ -1 -2 -2 0 -2]
[ 1 2 0 -2 2]
[ 2 1 -1 0 0]
sage: Y
[ 0 -2 0 1 -2]
[ 0 0 0 0 -1]
[ -1 -1 0 -1 2]
[ 0 -2 -2 -2 1]
[ -1 -1 1/2 0 0]
sage:
sage: parent(P)
Multivariate Polynomial Ring in x, y over Rational Field
sage:
sage: # do the substitution
sage: q = P.subs(x=X, y=Y)
sage: q
[-170339/256 221009/128 337381/128 90749/64 -27463/256]
[-131901/128 461983/64 574233/64 226241/32 -367357/128]
[-339165/128 441951/64 663751/64 184439/32 -89985/128]
[ 167421/128 -917935/64 -1076103/64 -456607/32 823585/128]
[ 8761/64 -3571/32 -6061/32 -815/8 3597/64]
sage: parent(q)
Full MatrixSpace of 5 by 5 dense matrices over Rational Field
sage: parent(q) is M
True
sage:
sage: # sanity check that it's what we expected
sage: q2 = sage_eval(str(P).replace('x','X').replace('y','Y'),locals=locals())
sage: q == q2
True