Hello, we do not have too much in the question. So here is my choice of an elliptic curve $E$ defined over my choice of a finite field $F$ of characteristic $p$. Of course, $E(F)$ is a finite group. For a small $n$ we can simply ask for the torsion points with that order. For bigger $n$ we may want to compute first the order of $E(F)$... This works only for small orders. Illustrating code/dialog with interrupting comments:

```
sage: p = 2020.next_prime()
sage: p
2027
sage: F = GF(p)
sage: E = EllipticCurve( F, [3, 4] )
sage: E.order()
2040
sage: E.order().factor()
2^3 * 3 * 5 * 17
```

so we expect no torsion point of order $7$.
For a "big $n$" we cannot use division polynomials.
But we could do this for $n$ among $2,3,4,5,6,8,10,15$ (even if $E$ has big order). Let us exemplify this.

For $n=2$ we can first ask for the corresponding division polynomial.

```
sage: E.torsion_polynomial(2)
4*x^3 + 12*x + 16
sage: E.torsion_polynomial(2).roots(multiplicities=0)
[2026, 1673, 355]
sage: [ E.lift_x(a).xy() for a in _ ] # _ is the above list
[(2026, 0), (1673, 0), (355, 0)]
```

We obtain as expected three $2$-torsion points.
(Well, an option `all=True`

is usually needed, we will see it later, but here...)
They correspond to the roots of the defining polynomial in $x$...

```
sage: R.<x> = PolynomialRing(F)
sage: f = x^3 + 3*x + 4
sage: f.factor()
(x + 1) * (x + 354) * (x + 1672)
sage: f.roots(multiplicities=False)
[2026, 1673, 355]
```

We have the order of $E(F)$, but we can also have its structure group:

```
sage: A = E.abelian_group()
sage: A
Additive abelian group isomorphic to Z/1020 + Z/2 embedded in Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 3*x + 4 over Finite Field of size 2027
```

So $E(F)\cong \Bbb Z/1020\oplus \Bbb Z/2$. Here, we expect $4$-torsion elements coming only from the $\Bbb/1020$-part, let us see them:

```
sage: G.<a,b> = AbelianGroup( [1020, 2] )
sage: G
Multiplicative Abelian group isomorphic to C1020 x C2
sage: [ g for g in G if g.order() == 4 ]
[a^255, a^255*b, a^765, a^765*b]
```

Here, we see four $4$-torsion elements, the $C_2\cong\Bbb Z/2$-component is either trivial or $b$, i.e. both chances are taken, the first component corresponds to $255$ (so thet $4\cdot 255=1020=0$ in that component) and $3\cdot 255$. We also expect four $4$-torsion elements in $E(F)$. Let us get them via the division polynomial.

```
sage: f = E.division_polynomial(4)
sage: f_roots = f.roots(multiplicities=False)
sage: points = []
sage: for a in f_roots:
....: try:
....: points.extend( E.lift_x(a, all=True) )
....: except Exception:
....: pass
....:
sage: [ P.xy() for P in points if P.order() == 4 ]
[(1978, 318), (1978, 1709), (759, 1046), (759, 981)]
```

We can obtain these points also by using the generators of $E$. Asking for and getting them is possible because of the reduced size of the data...

```
sage: P, Q = E.gens()
sage: P.order(), Q.order()
(1020, 2)
sage: [ (k*P + n*Q).xy() for k, n in cartesian_product( [[255, 765], [0,1]] ) ]
[(1978, 318), (759, 1046), (1978, 1709), (759, 981)]
```

Having $P,Q$ as above we can immediately get the points of order $3$, note that $1020/3 = 340$:

```
sage: [ (k*P).xy() for k in (340, 2*340) ]
[(633, 273), (633, 1754)]
```

The brutal way to get these points is by listing all of them and getting those of order exactly $3$:

```
sage: [ P.xy() for P in E.points() if P.order() == 3 ]
[(633, 273), (633, 1754)]
```

This works because of the reduced size of data.

To have a more pointed answer, please edit the question to have it with pointed details for the application in mind.