Numerics for elliptic Feynman integrals

The Standard Model involves several heavy particles: the Z- and W-bosons, the Higgs boson and the top quark. Precision studies of these particles require on the theoretical side quantum corrections at the two-loop order and beyond. It is a well-known fact that starting from two-loops Feynman integrals with massive particles can no longer be expressed in terms of multiple polylogarithms. This raises immediately the following question. What is the larger class of functions needed to express the relevant Feynman integrals? For single-scale two-loop Feynman integrals related to a single elliptic curve we have by now the answer: They are expressed as iterated integrals of modular forms [1]. This brings us to a second question: Is there an efficient method to evaluate these functions numerically in the full kinematic range? In this contribution we review how this can be done.