I can generate a random nxn matrix (for property based testing) by just by filling a nxn matrix with random numbers. But can anyone think of a way to generate random pairs of commutative matrices?

One way (but I’m not sure whether this is good or bad) is generate two random matrices A and B, and return (A+B,A-B). These two matrices certainly commute, as (A+B)(A-B) = (A-B)(A+B) = A^2 - B^2.