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(x) = \prod_{z : \varphi(z) = P}(x - 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:).