I recently implemented an algorithm to determine these number fields by computing the $j$-polynomial: Let $\varphi: X_0(N) \to E$ be a fixed modular parametrization and $P \in E(\mathbb{Q})$. By j-polynomial I mean the polynomial $F_P(j) = \prod_{z : \varphi(z) = P}(j - j(z))$. There's a Laurent series $x(q)$ with integer coefficients, which is the modular function $\frac{1}{x \circ \varphi}$ on $X_0(N)$,
We compute $x(q)$ via the gp function ell.taniyama. Then set $u = \frac{1}{j(q)}$, which is also an element in $\mathbb{Z}[[q]]$. Then using Linear algebra, one can find an irreducible polynomial $F$ such that $F(x,u) = 0$.
Setting $x_0 = \frac{1}{x(P)}$ and the polynomial $j^{2 \deg \varphi}F(x_0,1/j)$ is a constant multiple of $F_P(j)F_{-P}(j)$. Then we could use complex analytic method to determine which factor corresponds to $P$ and which is $-P$.
As an example we take the elliptic curve 121b1 with rank 1 and trivial torsion. $E(\mathbb{Q})$ is generated by $P = (4:5:1)$. Then we compute some j-polynomials:
$F_{-P}(x) = x^4 + 1421551441067913615636000 x^3 + 910640170936002098476853963114167004130307406250 x^2 - 55869041153225091624766256009488963324954953937500000 x + 1513370207928838475604980619812428055721351700525634765625$
$F_{4P}(x) = x^{4} - \frac{1131444376477476487694208}{43} x^{3} + \frac{11389877706351841520907948036498862509059802748293120}{1849} x^{2} + \frac{831545351967828972021160038394755202358001953700040409088}{1849} x + \frac{16095967144279358005293903881120972455827496828529236714192896}{1849}$
$F_{P}(x) = x^4 + 98823634118413525094400000 x^3 + 45688143672322270430861721600000000 x^2 - 496864268553728774541064273920000000000 x + 1577314437358442913340940353536000000000000$
Since atkin-lehner acts as +1 in this case the number fields in question are splitting fields of these polynomials, respectively.
I plan to compute with the rank 2 curve 389a mentioned in the original post. Doing so is hard with the current computing power I have, so I was thinking about replacing $u$ by an $\eta$-product with smaller valence.
Please let me know if this helps. I'll keep polishing this algorithm for the goal of including it in my coming up thesis:).