In this paper, we present an efficient algorithm to sample random sparse matrices to be used as check matrices for quantum Low-Density Parity-Check (LDPC) codes. To ease the treatment, we mainly describe our algorithm as a technique to sample a dual-containing binary LDPC code, hence, a sparse matrix $\mathbf H\in\mathbb F_2^{r\times n}$ such that $\mathbf H\mathbf H^\top = \mathbf 0$. However, as we show, the algorithm can be easily generalized to sample dual-containing LDPC codes over non binary finite fields as well as more general quantum stabilizer LDPC codes. While several constructions already exist, all of them are somewhat algebraic as they impose some specific property (e.g., the matrix being quasi-cyclic). Instead, our algorithm is purely combinatorial as we do not require anything apart from the rows of $\mathbf H$ being sparse enough. In this sense, we can think of our algorithm as a way to sample sparse, self-orthogonal matrices that are as random as possible. Our algorithm is conceptually very simple and, as a key ingredient, uses Information Set Decoding (ISD) to sample the rows of $\mathbf H$, one at a time. The use of ISD is fundamental as, without it, efficient sampling would not be feasible. We give a theoretical characterization of our algorithm, determining which ranges of parameters can be sampled as well as the expected computational complexity. Numerical simulations and benchmarks confirm the feasibility and efficiency of our approach.