Summary: Loop quantum cosmology is a symmetry-reduced application of loop quantum gravity that has led to the resolution of classical singularities such as the big bang, and those at the center of black holes. This can be seen through numerical simulations involving the quantum Hamiltonian constraint that is a partial difference equation. The equation allows one to study the evolution of sharply-peaked Gaussian wave packets that generically exhibit a quantum ‘bounce’ or a non-singular passage through the classical singularity, thus offering complete singularity resolution. In addition, von-Neumann stability analysis of the difference equation – treated as a stencil for a numerical solution that steps through the triad variables – yields useful constraints on the model and the allowed space of states. In this paper, we develop a new method for the numerical solution of loop quantum cosmology models using a set of basis functions that offer a number of advantages over computing a solution by stepping through the triad variables. We use the Corichi and Singh model for the Schwarzschild interior as the main case study in this effort. The main advantage of this new method is computational efficiency and the ease of parallelization. In addition, we also discuss how the stability analysis appears in the context of this new approach.