I will begin by introducing the Gaussian, Laguerre and Jacobi ensembles and their corresponding eigenvalue densities. The moments of these eigenvalue densities are generated by the corresponding resolvent, R(x). When investigating large matrices of size N, it is natural to expand R(x) as a series in 1/N, as N tends to infinity. The loop equation formalism enables one to compute R(x) to any desired order in 1/N via a triangular recursive system. This formalism is equivalent to the topological recursion, the Schwinger-Dyson equations and the Virasoro constraints, among other things. The loop equations provide a relatively accessible entry-point to these topics and my derivation will rely on nothing more than integration by parts, as Aomoto applied to the Selberg integral. Time permitting, I may also explore links to the topological recursion and/or some combinatorics. All original results will be from joint work with Peter Forrester and Nicholas Witte.