If we know the above answer, or generally if we know the field of definition for $E [ n ] (\overline{\mathbb Q})$, then we can initialize the curve $E$ in sage over this right field and ask for the corresponding $n^2$ torsion points of order dividing $n$. In the above case we have a simple situation:
sage: R.<x> = PolynomialRing(QQ)
sage: K.<J> = NumberField( x^2+x+1 )
sage: J^3
1
sage: L.<a> = K.extension( x^3 - 2 )
sage: a^3
2
sage: E = EllipticCurve( L, [0, -2] )
sage: Origin = E( 0 )
sage: Origin.division_points( 2 )
[(0 : 1 : 0), (a : 0 : 1), ((-J - 1)*a : 0 : 1), (J*a : 0 : 1)]
Note: It would have been possible to start with the simpler curve $y^2 = x^3 -1$, which is a twist of the given one.
Things get very complicated for bigger values of $n$.
For instance, if we want to get the points of order $3$ on the given curve, the extension needed is related to the multiplication by $3$ morphism on the curve. (I know, this is an instance of an other question. But the general path remains.) We may ask for it:
sage: E = EllipticCurve( QQ, [0,-2] )
sage: f,g = E.multiplication_by_m(3)
sage: f.factor()
(1/9) * x^-2 * (x - 2)^-2 * (x^2 + 2*x + 4)^-2 * (x^3 + 6*x^2 - 8) * (x^6 - 6*x^5 + 36*x^4 - 16*x^3 + 48*x^2 + 64)
sage: g.factor()
(1/27) * x^-3 * (x - 2)^-3 * y * (x^2 + 2*x + 4)^-3 * (x^3 + 16) * (x^3 - 6*x^2 - 12*x - 8) * (x^6 + 6*x^5 + 48*x^4 - 88*x^3 + 96*x^2 - 96*x + 64)
And a possibility to proceed would be to use them.
In fact, sage even offers us the relevant field.
EQ = EllipticCurve( QQ, [0, -2] )
K.<b> = EQ.division_field( 3 )
print K
EK = EQ.base_extend( K )
for P in EK.torsion_points():
print "Order %s :: %s" % ( P.order(), P )
And we get:
Order 3 :: (2 : 1/18*b^3 - 1/6*b^2 - 2/3*b + 7/9 : 1)
Order 3 :: (2 : -1/18*b^3 + 1/6*b^2 + 2/3*b - 7/9 : 1)
Order 3 :: (0 : 1/6*b^2 - 1/3*b - 1/3 : 1)
Order 3 :: (1/18*b^3 - 1/6*b^2 + 1/3*b - 11/9 : 1/18*b^3 - 1/6*b^2 - 2/3*b + 7/9 : 1)
Order 3 :: (-1/18*b^3 + 1/6*b^2 - 1/3*b - 7/9 : -1/18*b^3 + 1/6*b^2 + 2/3*b - 7/9 : 1)
Order 3 :: (0 : -1/6*b^2 + 1/3*b + 1/3 : 1)
Order 3 :: (-1/18*b^3 + 1/6*b^2 - 1/3*b - 7/9 : 1/18*b^3 - 1/6*b^2 - 2/3*b + 7/9 : 1)
Order 3 :: (1/18*b^3 - 1/6*b^2 + 1/3*b - 11/9 : -1/18*b^3 + 1/6*b^2 + 2/3*b - 7/9 : 1)
Order 1 :: (0 : 1 : 0)
(For the posted question one can replace the 3 in EQ.division_field( 3 ) by 2.)
The other possibility is of a numerical taste.
We may use the bijection from $\mathbb{C}/\Lambda$ to $E(\mathbb{C})$, usual notations, and compute to a good precision the image of the obvious $3$-torsion points (i.e. $\frac 13\cdot\Lambda$) in the domain of definition.
