Computing in this case means really computing. To see which is the complexity of the needed calculus, i am constrained to use some software helping me to get the answer quickly. When the answer stays here explicitly we can still argue if "easy / handy computations" may suffice.
In sage i am intializing the curve, asking for the $24$-division polynomial $f$ associated to it. In case of a nice information so far - we may still ask for the field where its roots live in, then make a base change and compute the generators (instead of computing all $24^2$ points). Well, instead of plotting the polynomial $f$ - which will soon have degree $289$, let us print the factors of $f$, or at least some information on them when they make the line explode to the right.
E = EllipticCurve(QQ, [-15, 22])
f = E.division_polynomial(24)
for fact, mul in f24.factor():
info = ( str(fact) if fact.degree() <= 8
else f'factor of degree {fact.degree()}' )
print(f'{info} with discriminant {ZZ(fact.discriminant()).factor()}')
We get so far:
x - 3 with discriminant 1
x - 2 with discriminant 1
x + 1 with discriminant 1
x^2 - 10*x + 13 with discriminant 2^4 * 3
x^2 - 4*x + 7 with discriminant -1 * 2^2 * 3
x^2 + 2*x - 11 with discriminant 2^4 * 3
x^3 - 9*x^2 + 51*x - 83 with discriminant -1 * 2^8 * 3^5
x^3 + 3*x^2 - 21*x + 25 with discriminant -1 * 2^8 * 3^3
x^4 + 4*x^3 - 66*x^2 + 148*x - 71 with discriminant 2^20 * 3^6
x^4 + 4*x^3 + 30*x^2 - 236*x + 313 with discriminant 2^20 * 3^7
x^6 + 18*x^5 - 141*x^4 + 428*x^3 - 1161*x^2 + 2658*x - 2507 with discriminant 2^44 * 3^19
x^8 - 40*x^7 + 268*x^6 - 952*x^5 + 1270*x^4 + 7592*x^3 - 40724*x^2 + 73400*x - 46703 with discriminant 2^88 * 3^30
x^8 - 16*x^7 + 172*x^6 - 592*x^5 + 454*x^4 - 2032*x^3 + 16876*x^2 - 37168*x + 25633 with discriminant 2^80 * 3^28
factor of degree 12 with discriminant 2^188 * 3^75
factor of degree 16 with discriminant 2^352 * 3^120
factor of degree 16 with discriminant 2^352 * 3^124
factor of degree 24 with discriminant 2^776 * 3^294
factor of degree 32 with discriminant 2^1408 * 3^504
factor of degree 48 with discriminant 2^3104 * 3^1164
factor of degree 96 with discriminant 2^12416 * 3^4632
And you really do not want to see (here) the last polynomial of degree $96$. One can compute of course closer information, so for the last question regarding a possible soft aimed to help, well sage is a free choice.
Here is as supplementary information what we can immediately compute, and may help for the asked purpose.
We can ask for the torsion points over $\Bbb Q$.
sage: [ T.xy() if T else 'O' for T in E.torsion_points() ]
[(-1, -6), (-1, 6), 'O', (2, 0), (3, -2), (3, 2)]
We have some polynomials of degree two listed as factors, discriminant is sometimes $3$ times a square, so it may be interesting to switch to the field $F=\Bbb Q(\sqrt3)=\Bbb Q(a)$, $a=\sqrt 3$. Then the given curve, seen over $F$ has...
R.<x> = PolynomialRing(QQ)
F.<a> = NumberField(x^2 - 3)
print(f'F is the field:\n{F}\n')
E = EllipticCurve(F, [-15, 22])
print(f'E has over F the torsion group:\n{E.torsion_subgroup()}\n')
print(f'Its generators are {E.torsion_subgroup().gens()}')
And the prints give the information:
F is the field:
Number Field in a with defining polynomial x^2 - 3
E has over F the torsion group:
Torsion Subgroup isomorphic to Z/6 + Z/2
associated to the Elliptic Curve defined by y^2 = x^3 + (-15)*x + 22
over Number Field in a with defining polynomial x^2 - 3
Its generators are ((-2a + 5 : -6a + 12 : 1), (2 : 0 : 1))
Here, instead of getting $E[24]$, let us compute starting from $F$ the pieces $E[3]\subset E(\overline{\Bbb Q})$ and $E[2^k]\subset E(\overline{\Bbb Q})$, $k=1,2,3$ as far as possible.
sage: E.division_polynomial(3).factor()
(3) * (x - 3) * (x^3 + 3*x^2 - 21*x + 25)
sage: E.division_polynomial(2).factor()
(4) * (x - 2a + 1) (x - 2) * (x + 2*a + 1)
sage: E.division_polynomial(4).factor()
(8) * (x - 2a + 1) (x - 2) * (x + 2a + 1) (x^2 + (-4a + 2)x + 8a - 11)
(x^2 - 4x + 7) (x^2 + (4a + 2)x - 8*a - 11)
and for the $8$-division let us better take each factor one by one...
sage: for fact, mul in E.division_polynomial(8).factor():
....: print(fact, 'with discriminant', F(fact.discriminant()).factor())
....:
x - 2*a + 1 with discriminant 1
x - 2 with discriminant 1
x + 2*a + 1 with discriminant 1
x^2 + (-4*a + 2)*x + 8*a - 11 with discriminant (15*a + 26) * (-a)^2 * (-a + 1)^8
x^2 - 4*x + 7 with discriminant (-4*a - 7) * (-a)^2 * (-a + 1)^4
x^2 + (4*a + 2)*x - 8*a - 11 with discriminant (209*a + 362) * (-a)^2 * (-a + 1)^8
x^4 + (-4*a - 8)*x^3 + (48*a + 78)*x^2 + (-132*a - 248)*x + 104*a + 241
with discriminant (296011017105*a + 512706121226) * (-a)^12 * (-a + 1)^36
x^4 + (4*a - 8)*x^3 + (-48*a + 78)*x^2 + (132*a - 248)*x - 104*a + 241
with discriminant (109552575*a + 189750626) * (-a)^12 * (-a + 1)^36
x^8 + (-16*a + 8)*x^7 + (128*a - 116)*x^6 + (-240*a - 232)*x^5
+ (-1376*a + 5398)*x^4 + (9232*a - 20104)*x^3 + (-24000*a + 31276)*x^2
+ (31472*a - 21016)*x - 17248*a + 5041
with discriminant
(319156261898201133032451783177661812885799840920*a
+ 552794861161443374798157973077388454680145589601)
* (-a)^56 * (-a + 1)^176
x^8 + (16*a + 8)*x^7 + (-128*a - 116)*x^6 + (240*a - 232)*x^5
+ (1376*a + 5398)*x^4 + (-9232*a - 20104)*x^3 + (24000*a + 31276)*x^2
+ (-31472*a - 21016)*x + 17248*a + 5041
with discriminant
(12011126751796371544078833423566810504717197292016904*a
+ 20803881790261051311370607096723572357290287974786657)
* (-a)^56 * (-a + 1)^176
sage:
Results were manually rearranged. Now let us exhibit the $3$-torsion. This is "hard enough" (to type). We ask for the corresponding division polynomial, which is of degree $3$, use $b$ for a "fixed symbolic root", so we expect to get a torsion point $T=(b,?)$, but the question mark involves taking the square root from $b^3 -15b+22$, denote it by $c$, so we have to extend $L=\Bbb Q(b)$ to $K=\Bbb Q(c)$. Over $K$ we have then at least a subgroup with $3^2$ elements with torsion order dividing $3$. Therin $3^2-1=8$ elements have order exactly $3$. Let us plot them in a final step...
Here is my dialog with the sage interpreter:
sage: E = EllipticCurve(QQ, [-15, 22])
sage: f = E.division_polynomial(3).factor()
sage: f
(3) * (x - 3) * (x^3 + 3*x^2 - 21*x + 25)
sage: L.<b> = NumberField(x^3 + 3x^2 - 21x + 25)
sage: (b^3 - 15*b + 22).is_square()
False
sage: # so we do not have an easy torsion point T = (b, ...) over L
sage: RL.<y> = PolynomialRing(L)
sage: K.<c> = L.extension(y^2 - (b^3 - 15*b + 22))
sage: EK = E.base_extend(K)
sage: EK.torsion_subgroup()
Torsion Subgroup isomorphic to Z/6 + Z/3 associated to the
Elliptic Curve defined by y^2 = x^3 + (-15)x + 22 over
Number Field in c with defining polynomial y^2 + 3b^2 - 6*b + 3
over its base field
sage: for T in EK.torsion_points():
....: if T.order() == 3:
....: print(T.xy())
....:
(3, 2)
(3, -2)
((3/8b^2 + 2b - 27/8)c - 1/2b - 3/2, (3/8b^2 + 3/2b - 55/8)c - 3/2b^2 - 15/2b + 18)
(b, c)
((-3/8b^2 - 2b + 27/8)c - 1/2b - 3/2, (3/8b^2 + 3/2b - 55/8)c + 3/2b^2 + 15/2b - 18)
((3/8b^2 + 2b - 27/8)c - 1/2b - 3/2, (-3/8b^2 - 3/2b + 55/8)c + 3/2b^2 + 15/2b - 18)
((-3/8b^2 - 2b + 27/8)c - 1/2b - 3/2, (-3/8b^2 - 3/2b + 55/8)c - 3/2b^2 - 15/2b + 18)
(b, -c)
In the last loop with lots of prints we see the two expected torsion points $(b,\pm c)$.