We study numerically the effect of sequence heterogeneity on the thermodynamic properties of a Poland-Scheraga model for DNA denaturation taking into account self-avoidance, i.e. with exponent c_p=2.15 for the loop length probability distribution. In complement to previous on-lattice Monte Carlo like studies, we consider here off-lattice numerical calculations for large sequence lengths, relying on efficient algorithmic methods. We investigate finite size effects with the definition of an appropriate intrinsic length scale x, depending on the parameters of the model. Based on the occurrence of large enough rare regions, for a given sequence length N, this study provides a qualitative picture for the finite size behavior, suggesting that the effect of disorder could be sensed only with sequence lengths diverging exponentially with x. We further look in detail at average quantities for the particular case x=1.3, ensuring through this parameter choice the correspondence between the off-lattice and the on-lattice studies. Taken together, the various results can be cast in a coherent picture with a crossover between a nearly pure system like behavior for small sizes N < 1000, as observed in the on-lattice simulations, and the apparent asymptotic behavior indicative of disorder relevance, with an (average) correlation length exponent \nu_r >= 2/d (=2).