  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.
