We present a numerical method for the evaluation of the mass gap, and the low-lying energy gaps, of a large family of free-fermionic and free-parafermionic quantum chains. The method is suitable for some generalizations of the quantum Ising and XY models with multispin interactions. We illustrate the method by considering the Ising quantum chains with uniform and random coupling constants. The mass gaps of these quantum chains are obtained from the largest root of a characteristic polynomial. We also show that the Laguerre bound, for the largest root of a polynomial, used as an initial guess for the largest root in the method, gives us estimates for the mass gaps sharing the same leading finite-size behavior as the exact results. This opens an interesting possibility of obtaining precise critical properties very efficiently which we explore by studying the critical point and the paramagnetic Griffiths phase of the quantum Ising chain with random couplings. In this last phase, we obtain the effective dynamical critical exponent as a function of the distance-to-criticality. Finally, we compare the mass gap estimates derived from the Laguerre bound and the strong-disorder renormalization-group method. Both estimates require comparable computational efforts, with the former having the advantage of being more accurate and also being applicable away from infinite-randomness fixed points. We believe this method is a new and relevant tool for tackling critical quantum chains with and without quenched disorder.