  We numerically study the binary disordered Poland-Scheraga model of DNA
denaturation, in the regime where the pure model displays a first order
transition (loop exponent $c=2.15>2$). We use a Fixman-Freire scheme for the
entropy of loops and consider chain length up to $N=4 \cdot 10^5$, with
averages over $10^4$ samples. We present in parallel the results of various
observables for two boundary conditions, namely bound-bound (bb) and
bound-unbound (bu), because they present very different finite-size behaviors,
both in the pure case and in the disordered case. Our main conclusion is that
the transition remains first order in the disordered case: in the (bu) case,
the disorder averaged energy and contact densities present crossings for
different values of $N$ without rescaling. In addition, we obtain that these
disorder averaged observables do not satisfy finite size scaling, as a
consequence of strong sample to sample fluctuations of the pseudo-critical
temperature. For a given sample, we propose a procedure to identify its
pseudo-critical temperature, and show that this sample then obeys first order
transition finite size scaling behavior. Finally, we obtain that the disorder
averaged critical loop distribution is still governed by $P(l) \sim 1/l^c$ in
the regime $l \ll N$, as in the pure case.
