Abstract: Building on the work of Iwaniec, Luo and Sarnak, we use the $n$-level density to bound the probability of vanishing to order at least $r$ at the central point for families of cuspidal newforms of prime level $N \to \infty$, split by sign. There are three methods to improve bounds on the order of vanishing: optimizing the test functions, increasing the support, and increasing the $n$-level density studied. Previous work determined the optimal test functions for the $1$ and $2$-level densities in certain support ranges, leading to marginal improvements in bounds and making it not a productive avenue for further research. Similarly the support has been increased as far as possible, and further progress is shown to be related to delicate and difficult conjectures in number theory. Thus we concentrate on the third method, and study the higher centered moments (which are similar to the $n$-level densities but combinatorially easier). We find the level at each rank for which the bounds on the order of vanishing is the best, thus producing world-record bounds on the order of vanishing to rank at least $r$ for every $r > 2$ (for example, our bounds for vanishing to order at least 5 or at least 6 are less than half the previous bounds, a significant improvement). Additionally, we explicitly calculate the optimal test function for the $1$-level density from previous work and compare it to the naive test functions for higher levels. In doing so, we find that the optimal test function for certain levels are not the optimal for other levels, and some test functions may outperform others for some levels but not in others. Finally, we explicitly calculate the integrals needed to determine the bounds, doing so by transforming an $n$-dimensional integral to a $1$-dimensional integral and greatly reducing the computation cost in the process.