The modularity theorem (original Shimura-Taniyama-Weil conjecture) asserts the existence of a covering (uniformization) map $\pi:X_0(N) \to E$ for every $E$, an elliptic curve defined over $\mathbb{Q}$.
If $y^2 = x^3 + a x + b$ is the Weierstraß equation of $E$, this amounts to finding power-series $x(q)$, $y(q)$ in $q = \exp(2 \pi i z)$ where $z$ is the standard coordinate of $\mathbb{C}$. The $x(q), y(q)$ are from the function field of $X_0(N)$, so they are modular functions with respect to $\Gamma_0(N)$ (is this true?).
In the system Pari/GP there is the function elltaniyama to compute such $x(q)$, $y(q)$, and in Sage there is the function modular_parametrization. But all searching with Google did not find me a paper, where the algorithm involved is described. I found some articles describing how to compute the modularity conductor $N$, but not, how $x(q)$ and $y(q)$ can be calculated.
Could someone direct me to a paper (or book) where to find an explanation of the algorithm used?
I don't know the inner workings of those computer algebra functions, but here's how you can explicitly compute the modular parametrization $\newcommand{\C}{\mathbb C} \renewcommand{\H}{\mathfrak H} \varphi: X_0(N) \to E$.
The unique (up to scaling) holomorphic differential $\omega$ on $E$ pulls back under $\varphi$ to a differential $\varphi^*\omega$ on $X_0(N)$. Thinking of $X_0(N)$ as $\Gamma_0(N) \backslash \H^*$, we can consider $\varphi^*\omega$ as a $\Gamma_0(N)$-invariant differential on $\H^*$. This can be written as $$ \varphi^* \omega = 2\pi i f(\tau) \, d\tau $$ where $f$ is a weight $2$ cusp form. (Indeed, the one corresponding to $E$ in other equivalent statements of modularity.) Letting $q = e^{2\pi i \tau}$ and writing $f(\tau) = \sum_{n=1}^\infty a_n(f) q^n$ in its $q$-expansion, we can integrate the differential $$ 2\pi i f(\tau) \, d\tau = f(q) \frac{dq}{q} $$ to get a map \begin{align*} \H^* &\to \C/\Lambda\\ \tau_0 &\mapsto 2\pi i \int_{\tau_0}^\infty f(\tau) \, d\tau = \int_0^{q_0} f(q) \frac{dq}{q} \end{align*} where $q_0 = e^{2\pi i \tau_0}$ and $\Lambda \leq \C$ is the period lattice of $E$. To write the elliptic curve $\C/\Lambda$ in Weierstrass form, we apply the inverse Abel-Jacobi map (sometimes called the elliptic exponential) using the Weierstrass $\wp$-function. This is the map \begin{align*} \C/\Lambda &\to E: y^2 = x^3 + a x + b\\ z &\mapsto (x,y) = \left(\wp(z), \frac{1}{2} \wp'(z)\right) \, , \end{align*} and composing the two maps give the series $x(q), y(q)$.
A few notes on computational effectivity. The Fourier coefficients $a_n(f)$ can be computed by counting points. Assuming $E$ is in minimal Weierstrass form, then $\#E(\mathbb{F}_p) = p + 1 - a_p(f)$ for all primes $p \nmid N$. (Here we are again using modularity, in that $a_p(f) = a_p(E)$.) There are efficient methods for counting points, such as Schoof's algorithm and its variants. Using the formulas \begin{align*} a_{mn} &= a_m a_n \text{ if $\gcd(m,n)=1$}\\ a_{p^n} &= a_p a_{p^{n-1}} - p a_{p^{n-2}} \end{align*} we can recover all the Fourier coefficients $a_n = a_n(f)$.
The values of $\wp(z)$ can also be computed effectively. Write $\wp(z) = \frac{1}{z^2} + \frac{c_{-1}}{z} + \sum_{n=0}^\infty c_n z^n$ as a Laurent series. Similarly writing $\frac{1}{2} \wp'(z)$ as a Laurent series and substituting in to the equation for $E$ \begin{align*} \left(\frac{\wp'(z)}{2}\right)^2 = (\wp(z))^3 + a\wp(z) + b \end{align*} allows us to find a recursive formula for the coefficients $c_n$. (The actual formula is pretty hideous, but it can be done.)
For example, consider the elliptic curve $E: y^2 = x^3 + 2x - 3$, curve 440.b2 in the LMFDB. As shown on that page, the corresponding modular form $f$ has $q$-expansion $$ f(q) = q - q^{5} - 2q^{7} - 3q^{9} + q^{11} - 4q^{13} - 4q^{17} + O(q^{20}) \,. $$ Then $$ F(q) := \int f(q) \frac{dq}{q} = q - \frac{1}{5}q^{5} - \frac{2}{7}q^{7} - \frac{1}{3}q^{9} + \frac{1}{11}q^{11} - \frac{4}{13}q^{13} - \frac{4}{17}q^{17} + O(q^{20}) \, . $$ Using the recurrence described above, we find the Laurent series \begin{align*} \wp(z) &= \frac{1}{z^{2}} - \frac{2}{5}z^{2} + \frac{3}{7}z^{4} + \frac{4}{75}z^{6} - \frac{18}{385}z^{8} + \frac{2591}{238875}z^{10} + \frac{8}{1925}z^{12}\\ &\qquad - \frac{72376}{44669625}z^{14} + \frac{159}{6131125}z^{16} + O(z^{17}) \end{align*} Composing these series and letting $x(q) = \wp(F(q))$ and $y(q) = \wp'(F(q))/2$ we find \begin{align*} x(q) &= q^{-2} + q^{4} + q^{6} + q^{10} + 2q^{14} + q^{16} + O(q^{17})\\ y(q) &= q^{-3} - q - q^{5} - q^{7} - 3q^{9} - 2q^{11} - 2q^{13} - 2q^{15} + O(q^{16}) \end{align*} which agrees with Sage: