y-coordinate of a 4-torsion point

edit

I would recommend to avoid Sage's symbolic ring for such computations.

The best here is to use algebraic numbers, which form the ring QQbar, and to work with polynomials over that ring.

For reference here is the Sage version I am working with.

sage: version()
'SageMath version 7.3, Release Date: 2016-08-04'

This is one way you could work.

sage: E3 = EllipticCurve([0, 0, 0, -3267, 45630])
sage: E = E3.division_polynomial(4); E
8*x^9 - 156816*x^7 + 7665840*x^6 - 25044299280*x^4 + 1873623790224*x^3 - 35065596749040*x^2 - 258340569679368*x + 6648422400893520
sage: E.parent()
Univariate Polynomial Ring in x over Rational Field

Here is how this polynomial factors over QQ:

sage: E.factor()
(8) * (x - 15) * (x^2 - 30*x + 2817) * (x^2 + 15*x - 3042) * (x^4 + 30*x^3 - 18252*x^2 + 280530*x + 6465339)

The list of factors over QQ:

sage: list(E.factor())
[(x - 15, 1),
 (x^2 - 30*x + 2817, 1),
 (x^2 + 15*x - 3042, 1),
 (x^4 + 30*x^3 - 18252*x^2 + 280530*x + 6465339, 1)]

Let us select one factor:

sage: p = E.factor()[1][0]; p
x^2 - 30*x + 2817
sage: p.parent()
Univariate Polynomial Ring in x over Rational Field

This factor lives in QQ[x]. Let us change to QQbar[x].

sage: P.<x> = QQbar[]
sage: p = P(p); p
x^2 - 30*x + 2817
sage: p.parent()
Univariate Polynomial Ring in x over Algebraic Field

Now this polynomial factors completely, and has two roots of multiplicity one.

sage: p.factor()
(x - 15.00000000000000? - 50.91168824543143?*I) * (x - 15.00000000000000? + 50.91168824543143?*I)
sage: p.roots()
[(15.00000000000000? - 50.91168824543143?*I, 1),
 (15.00000000000000? + 50.91168824543143?*I, 1)]

Forgetting about multiplicities:

sage: roots = p.roots(multiplicities=False); roots
[15.00000000000000? - 50.91168824543143?*I,
 15.00000000000000? + 50.91168824543143?*I]

Let us enter by hand the polynomial associated with E3.

sage: q = x^3 - 3267*x + 45630; q
x^3 - 3267*x + 45630
sage: q.parent()
Univariate Polynomial Ring in x over Algebraic Field

[Exercise: try to get hold of this polynomials starting from E3 instead of entering it by hand.]

Now for each root r of p, the corresponding y-value on E3 is sqrt(q(r)).

sage: for (x, y) in [(r, sqrt(q(r))) for r in roots]:
....:     print 'xr = {} = {}'.format(x, x.radical_expression())
....:     print 'yr = {} = {}\n'.format(y, y.radical_expression())
....:     
xr = 15.00000000000000? - 50.91168824543143?*I = -36*I*sqrt(2) + 15
yr = 293.1820459230292? + 450.1063341607327?*I = 27*sqrt(256*I*sqrt(2) - 160)

xr = 15.00000000000000? + 50.91168824543143?*I = 36*I*sqrt(2) + 15
yr = 293.1820459230292? - 450.1063341607327?*I = 27*sqrt(-256*I*sqrt(2) - 160)

These are algebraic numbers and maybe their minimal polynomial is useful.

sage: for (x, y) in [(r, sqrt(q(r))) for r in roots]:
....:     print 'xr has minpoly {}'.format(x.minpoly())
....:     print 'yr has minpoly {}\n'.format(y.minpoly())
....:     
xr has minpoly x^2 - 30*x + 2817
yr has minpoly x^4 + 233280*x ...