We present a parallel modular algorithm for finding characteristic polynomials of matrices with integer coefficient bivariate monomials. For each prime, evaluation and interpolation gives us the bridge between polynomial matrices and matrices over a finite field so that the Hessenberg algorithm can be used. After optimizations, we are able to save a significant amount of work by incremental Chinese remaindering and early termination.