As Stankewicz explains, although elliptic curves appear in two guises in the modular parameterization $X_0(N) \to E,$ first because $E$ is an elliptic curve, and secondly because $X_0(N)$ parameterizes elliptic curves, it is something of a red herring to think of these two appearances of elliptic curves as having anything to do with one another.
The reason that $X_0(N)$ appears in the problem of describing elliptic curves is because
elliptic curves have two dimensional $H^1$, and $X_0(N)$ is the Shimura variety over $\mathbb Q$ associated to the group $GL_2$. Thus, as Stankewicz notes, Shimura curves (which are the Shimura varieties attached to twisted forms of $GL_2$) can equally well give parameterizations of elliptic curves.
Now the way we prove things about $X_0(N)$ (e.g. properties of its Heegner points, as in Pete Clark's answer) is using its moduli interpretation. But there are two things to bear in mind:
First, most (maybe all?) properties of the special points on $X_0(N)$, such as the Heegner
points, are special cases of general aspects of the theory of Shimura varieties (so although the proofs use the moduli interpretation, the statements can be formulated in a way that
doesn't refer to the moduli-theoretic interpretation, but instead refers to the interpretation
of $X_0(N)$ as a Shimura variety).
Second, the transfer of information is always from $X_0(N)$ to $E$. So while Heegner points give certain interesting points on $X_0(N)$ defined over class fields of quadratic imaginary fields, which can be mapped down to $E$ to give interesting points on $E$ defined over such fields, if one takes a random point on $E$ defined over a class field of an imaginary quadratic field and pulls it back to $X_0(N)$, it is not so easy to say what is going on with the preimages in general.
Finally, I think remark (3) in Pete Clark's answer is an interesting one. In the Mazur and Swinnerton-Dyer paper that David Hansen refer's to in his first comment, if I am remembering correctly, they also suggest that the images in $E$ of the critical points of the map from $X_0(N)$ to $E$ that lie on the geodesic arc joining $0$ to $\infty$ in the upper half-plane may be worth studying. As with Birch's suggestion of Weierstrass points, I'm not sure how much has been done on this.
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:).
Best Answer
One expects that the majority of algebraic curves over number fields having genus $> 1$ should not be modular in this sense.
For instance, take a sufficiently general genus 2 curve $C$ over $\mathbf{Q}$. Then its $\ell$-adic $H^1$ (which is just the $\ell$-adic Tate module of its Jacobian) will be a 4-dimensional Galois representation whose image lands inside $\mathrm{GSp}_4(\mathbf{Z}_\ell)$. If $C$ is sufficiently generic, then the image of this Galois representation should be the whole of $\mathrm{GSp}_4(\mathbf{Z}_\ell)$ for all but finitely many $\ell$; in particular, it will be absolutely irreducible. (I don't know if this is known, but certainly one expects it to be the case.) On the other hand, if $C$ admits a non-constant map from $X_0(N)$, then the its $H^1$ would have to be a quotient of the $H^1$ of $X_0(N)$, and this can be calculated in terms of modular forms; in particular all its absolutely irreducible subquotients have dimension 2. So most genus 2 curves $C$ will not be modular in your sense, and if you get one that is, you should regard it as a rather unlikely coincidence.
(A more high-powered interpretation of this is that $H^1(C)$ should be the Galois representation attached to a degree 2 Siegel modular form. In some very special cases this Siegel modular form will be endoscopic, i.e. describable in terms of lifts from elliptic modular forms, but most Siegel mod forms will not be endoscopic and thus will not have anything to do with $X_0(N)$ for any $N$.)
If you're willing to relax your definition of "modular", though, you can get many more possibilities. There's a very striking result of Belyi stating that any algebraic curve defined over a number field can be obtained as the quotient of the upper half-plane by some subgroup of $PSL(2, \mathbf{Z})$, although the corresponding group will usually not be a congruence subgroup.