The analytical calculation of Feynman integrals is an important problem in modern quantum field theory. This task is important both for obtaining the most accurate predictions for observable quantities and for some areas of pure mathematics such as theory of periods in algebraic geometry. Nevertheless, it is not always possible to obtain analytical solutions for all Feynman integrals beyond one loop. All difficulties are usually associated with the appearance of elliptic or more complex geometric structures which inevitably arise when taking into account the masses of propagators.
In this work, we use an example of two loop elliptic master integrals arising from non-relativistic QCD as a laboratory to develop new methods for calculating non-polylogarithmic Feynman integrals. First of all, we will consider a new method that allows to obtain exact, in terms of the dimensional regularization parameter, solutions for the integrals under consideration. In this case, the solutions are expressed in terms of well-converging Frobenius power series. We will also briefly consider a new method for obtaining integral representations for the same integrals.