As per @Alexei's suggestion to use HoldForm
, here's a function that formats a polynomial the way I was looking for. In this case it's for hard-coded 3 variables, but could probably be straightforwardly extended to use the length of var
to do so automatically:
formatPolynomial[poly_, vars_ : {p, q, pdq}] := Module[{ dim, coeffs, output },
dim = Max[Cases[CoefficientRules[poly, vars], v_?VectorQ :> Total[v], 2]] + 1;
coeffs = CoefficientList[poly, vars, {dim, dim, dim}];
output = Sum[coeffs[[i, j, k]] Subscript[tmp, i - 1, j - 1, k - 1], {i, 1,
dim, 1}, {j, 1, dim, 1}, {k, 1, dim, 1}];
output = output /. {a_ Subscript[tmp, n_, m_, l_] :> HoldForm[a] Subscript[tmp, n, m,l]} /. {Subscript[tmp, n_, m_, l_] :> vars[[1]]^n vars[[2]]^m vars[[3]]^l} /. {a_ HoldForm[b_] :> HoldForm[({b} // MatrixForm)] a};
output
]
This gives output so the coefficients of individual monomials (which I'm looking for) are easily identifiable:

Apart
? $\endgroup$