We present a general approach to numerically compute the solutions of Riemann-Hilbert problems with jump conditions supported on disjoint intervals. Applied to the Fokas-Its-Kitaev Riemann-Hilbert problem, this enables the computation of orthogonal polynomials on multiple intervals, requiring only O(N) arithmetic operations to compute the first N recurrence coefficients. Such recurrence coefficients describe the flow of a semi-infinite Toda lattice, and expansions in these orthogonal polynomials yield a novel iterative method for solving indefinite linear systems and computing matrix functions, further yielding a fast algorithm for computing the solutions of Sylvester matrix equations. Other Riemann-Hilbert problems of this form yield finite-genus and soliton gas solutions of the Korteweg-de Vries equation. In particular, we compute large-genus solutions to simulate dispersive quantization and evaluate soliton gas solutions before asymptotic estimates are valid.