We calculate 3-loop master integrals for heavy quark correlators and the 3-loop quantum chromodynamics corrections to the ρ-parameter. They obey non-factorizing differential equations of second order with more than three singularities, which cannot be factorized in Mellin-N space either. The solution of the homogeneous equations is possible in terms of 2F1 Gauß hypergeometric functions at rational argument. In some cases, integrals of this type can be mapped to complete elliptic integrals at rational argument. This class of functions appears to be the next one arising in the calculation of more complicated Feynman integrals following the harmonic polylogarithms, generalized polylogarithms, cyclotomic harmonic polylogarithms, square-root valued iterated integrals, and combinations thereof, which appear in simpler cases. The inhomogeneous solution of the corresponding differential equations can be given in terms of iterative integrals, where the new innermost letter itself is not an iterative integral. A new class of iterative integrals is introduced containing letters in which (multiple) definite integrals appear as factors. For the elliptic case, we also derive the solution in terms of integrals over modular functions and also modular forms, using q-product and series representations implied by Jacobi’s ϑi functions and Dedekind’s η-function. The corresponding representations can be traced back to polynomials out of Lambert–Eisenstein series, having representations also as elliptic polylogarithms, a q-factorial 1/ηk(τ), logarithms, and polylogarithms of q and their q-integrals. Due to the specific form of the physical variable x(q) for different processes, different representations do usually appear. Numerical results are also presented.