I revisit the problem of solving multi-matrix systems through numerical large N methods. The framework is the constrained optimization of a loop truncated, collective, loop space representation of these systems. The loop space constraints are handled by the use of master variables. For the large N dynamics of two massless Yang-Mills coupled matrix quantum mechanics, the method is successfully applied directly in the massless limit for a range of values of the Yang- Mills coupling constant, and the scaling behaviour of different physical quantities derived from their dimensions are obtained with a high level of precision. We consider both planar properties of the theory, such as the large N ground state energy and multi-matrix correlator expectation values, and also the spectrum of the theory. For the spectrum, we establish that the U (N ) traced fundamental constituents remain massless and decoupled from other states, and that bound states develop well defined mass gaps, with the mass of the two degenerate lowest lying bound states being determined with a particularly high degree of accuracy.