In this paper, an analogue of matrix models from free probability is developed in the bi-free setting. A bi-matrix model is not simply a pair of matrix models, but a pair of matrix models where one element in the pair acts by left-multiplication on matrices and the other element acts via a `twisted'-right action. The asymptotic distributions of bi-matrix models of Gaussian random variables tend to bi-free central limit distributions with certain covariance matrices. Furthermore, many classical random matrix results immediately generalize to the bi-free setting. For example, bi-matrix models of left and right creation and annihilation operators on a Fock space have joint distributions equal to left and right creation and annihilation operators on a Fock space and are bi-freely independent from the left and right action of scalar matrices. Similar results hold for bi-matrix models of $q$-deformed left and right creation and annihilation operators provided asymptotic limits are considered. Finally bi-matrix models with asymptotic limits equal to Boolean independent central limit distributions and monotonically independent central limit distributions are constructed.