We calculate 3-loop master integrals for heavy quark correlators and the 3-loop QCD corrections to the rho-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 convergent close integer power series as 2_F_1 Gauss 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 theta_i functions and Dedekind's eta-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/eta^k(\tau), 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. The present results cover classes of elliptic solutions of a more general kind than those of the well-known sunrise and kite diagrams. Numerical results for the solutions are also presented.