Lyapunov functions provide a tool to analyze the stability of nonlinear systems without extensively solving the dynamics. Recent advances in sum-of-squares methods have enabled the algorithmic computation of Lyapunov functions for polynomial systems. However, for general large-scale nonlinear networks it is yet very difficult, and often impossible, both computationally and analytically, to find Lyapunov functions. In such cases, a system decomposition coupled to a vector Lyapunov functions approach provides a feasible alternative by analyzing the stability of the nonlinear network through a reduced-order comparison system. However, finding such a comparison system is not trivial and often, for a nonlinear network, there does not exist a single comparison system. In this work, we propose a multiple comparison systems approach for the algorithmic stability analysis of nonlinear systems. Using sum of-squares methods we design a scalable and distributed algorithm which enables the computation of comparison systems using only communications between the neighboring subsystems. We demonstrate the algorithm by applying it to an arbitrarily generated network of interacting Van der Pol oscillators. (c) 2016 Elsevier Ltd. All rights reserved.