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:
sage: E = EllipticCurve([0, 0, 0, 2, -3])
sage: x, y = E.modular_parametrization().power_series()
sage: x
q^-2 + q^4 + q^6 + q^10 + 2*q^14 + q^16 + O(q^18)
sage: y
-q^-3 - q - q^5 - q^7 - 3*q^9 - 2*q^11 - 2*q^13 - 2*q^15 + O(q^17)
