We describe an algorithm that we used to compute the q-expansions of all weight 2 cusp forms of prime level at most 2,000,000 and dimension at most 6. We also present an algorithm that we used to verify that there was only one cusp form of dimension 7 or more per Atkin-Lehner eigenspace for prime levels between 10,000 and 1,000,000. Our algorithm is based on Mestre's M\'ethode des Graphes, and involves supersingular isogeny graphs and Wiedemann's algorithm for finding the minimal polynomial of sparse matrices over finite fields.
Comment: v3: Updated with recent data