We develop an efficient numerical method to study the quantum critical behavior of disordered systems with O(N ) order-parameter symmetry in the large−N limit. It is based on the iterative solution of the large−N saddle-point equations combined with a fast algorithm for inverting the arising large sparse random matrices. As an example, we consider the superconductor-metal quantum phase transition in disordered nanowires. We study the behavior of various observables near the quantum phase transition. Our results agree with recent renormalization group predictions, i.e., the transition is governed by an infinite-randomness critical point, accompanied by quantum Griffiths singularities. In contrast to the existing numerical approach to this problem, our method gives direct access to the temperature dependencies of observables. Moreover, our algorithm is highly efficient because the numerical effort for each iteration scales linearly with the system size. This allows us to study larger systems, with up to 1024 sites, than previous methods. We also discuss generalizations to higher dimensions and other systems including the itinerant antiferromagnetic transitions in disordered metals.Copyright line will be provided by the publisher 1 Introduction Randomness can have much more dramatic effects at quantum phase transitions than at classical phase transitions because quenched disorder is perfectly correlated in the imaginary time direction which needs to be included at quantum phase transitions. Imaginary time acts as an additional coordinate with infinite extension at absolute zero temperature. Therefore, the impurities and defects are effectively very large which leads to strong-disorder phenomena including power-law quantum Griffiths singularities [1,2,3], infinite-randomness critical points characterized by exponential scaling [4,5], and smeared phase transitions [6]. For example, the zero-temperature quantum phase transition in the random transverse-field Ising model is governed by an infiniterandomness critical point [5] featuring slow activated (exponential) rather than power-law dynamical scaling. It is accompanied by quantum Griffiths singularities. This means, observables are expected to be singular not just at criticality but in a whole parameter region near the quan-