How do you explicitly compute the p-torsion points on a general elliptic curve in Weierstrass form? Consider the Weierstrass cubic
$$y^2z = x^3 + A\, xz^2+B\,z^3.$$
This defines a curve $E$ in $\mathbb{P}^2$, which if smooth is an elliptic curve with basepoint at $[0,1,0]$.
I'm interested in having an explicit description of the locus of $p$-torsion points of this curve, where $p$ is prime.  
In fact, suppose $p\neq 3$.  Then ideally I'd like to be able to find a curve $C$ in $\mathbb{P}^2$, given by an equation $f=0$ of degree $d=(p^2-1)/3$, so that the scheme $X=E\times_{\mathbb{P}^2} C$ is precisely the locus of points of exact order $p$.  
Example: For $p=2$, it's well known that $f=y$ gives such a curve.
I'd like $f$ to be an expression which depends on $A$ and $B$; i.e., I want to do this over a generic part of the moduli stack.  I would also like this expression to work in characteristic p;  in this case, $X$ should turn out to be the "scheme representing Drinfeld level structures $\mathbb{Z}/p\to E$".  (Edit: I'm particularly interested in families of curves which include supersingular curves.)
(My example curve $E$ is never smooth in characteristic $2$, but if you consider a more general Weierstrass form which is smooth in char. $2$, then you can find a degree $1$ curve $C$ which does what I ask.  For instance, if $E: y^2z+A\,xyz+yz^2=x^3$, take $f=A\,x+2\,y+z$.)
So my questions are:


*

*Is it usually possible to find an equation $f=0$ such that $E\cap C$ is exactly the $p$-torsion?  (Is this the same as asking that $X$ is a complete intersection?)  Can you ever show it's not possible?

*Are there known methods for computing the locus of $p$-torsion points explicitly?  Are there software packages which do this?  (I'm aware there are ways to find explicit torsion points on elliptic curves defined over some field or number ring; I'm asking for something a little different, I think.)

*Have people carried out these sorts of computations for various small values of $p$ (even $p=5$), and are these computations described in print?  (I'm probably most interested in this question.)  


Warning: I am not an algebraic geometer or number theorist.
 A: Since you asked about software, I'd just like to point out (if you don't know already) that SAGE (available at sagemath.org) can compute division polynomials easily.  The commands
R.<A,B> = PolynomialRing(GF(5)) 
E = EllipticCurve([A,B])
f = E.division_polynomial(5)
f
return the result
2*A*x^10 - A^2*B*x^5 + A^6 - 2*A^3*B^2 - B^4.
The warnings of the commenters apply;  you must be cautious about interpreting the division polynomials.  It is true that the roots of this polynomial give you the "physical" 5-torsion points of the elliptic curve in characteristic 5, but that's about all it says.  The polynomial does not tell you, eg, the structure of the group scheme E[5] over the supersingular locus.  
A little more helpful might be the formal group associated to this family:
G = E.formal_group()
G.mult_by_n(5,30)
which returns
2*A*t^5 + (2*A^6 - A^3*B^2 - B^4)*t^25 + O(t^30)
There's also a command G.group_law() whose output I'm not going to give here.  Together these data give you the structure of E[5] over an infinitesimal neighborhood of the zero section of your elliptic curve over $\text{Spec}\mathbf{F}_5[A,B,\Delta^{-1}]$. 
I'm a firm believer in explicit calculations as a means of developing intuition about algebro-geometric concepts.  But by all means read Katz-Mazur :-).
A: http://en.wikipedia.org/wiki/Division_polynomials
That's not the best wikipedia page. "The division polynomials form an elliptic divisibility sequence." is mentioned well before the far more important "the roots of the n'th division polynomial tell you the n-torsion in the curve".
A: Here is an attempt to answer my own question, using the "division
polynomials" of Kevin's and Jared's answers.   It is probably the maximally naive idea, and I do not claim it
works, though it's not clear to me that it can't.
I've community wiki-ed this answer, as it's probably junk anyway ...
Fix $A,B\in \mathbb{Z}$, obtaining an elliptic curve $E$ over $S=Spec
\mathbb{Z}[\Delta^{-1}]$, and a prime $p\geq 5$ not dividing the
discriminant $\Delta=\Delta(A,B)$.  The
division polynomial $\psi_p(t)$ is supposed to be a polynomial of
degree $d=(p^2-1)/2$ over $\mathbb{Z}$, whose roots are the values 
$t(P)=x(P)/z(P)$ as $P$ ranges over the points of exact order $p$ in
$E$.
Turn $\psi_p$ into a homogeneous polynomial $g$ of degree $d$ in
$\mathbb{Z}[x,y,z]$, so that $g(x,y,1)=\psi_p(x)$.  The polynomial
deterimes a curve $C=(g)$ in $P^2/S$, and thus a closed subscheme
$D=E\cap C$ of $E$.  Over $\mathbb{Z}[\Delta^{-1},p^{-1}]$, $D$ should
consist of 
the points of exact order $p$ on $E$ (with multiplicity $1$), together
with the basepoint of $E$ with multiplicity $d$.  
Claim. $D$ is an effective Cartier
divisor on $E/S$, of degree $3d$.  
Proof.  I don't know.  I need to prove that $D\to
S$ is flat, the main concern being the behavior over
$\mathbb{Z}_{(p)}$.  I don't even know if this is really plausible in general.  
Let's pretend we somehow know $D$ is an effective Cartier divisor on $E$
relative to the base $S$.  There is another relative Cartier
divisor, namely
$$ D' = E[p] \quad + \quad (d-1)[0]$$
where $E[p]$ is the $p$-torsion subgroup scheme of $E$, and $[0]$ is the
degree one divisor of the basepoint of $E$.
It seems clear that away from characteristic $p$, the divisors $D$
and $D'$ are equal.  Equality of effective divisors on a smooth curve
is a closed condition, so they should be equal over all of
$S$.  
The divisor I want is thus $D''=D-d[0]=D'-d[0]$ (which is still
effective).  Then it's really easy to find a homogeneous polynomial
$h$ of degree $2d/3$ which defines $D''$; because of the form of the
Weierstrass equation, you can produce it from $g$ by hand, and if my claim is true you can produce it globally, i.e., with coefficients in $\mathbb{Z}[A,B]$.
I've carried out this out in the case $p=5$, and it appears to "work".
That is, I get an answer which appears sane for general $A$ and $B$, and
which for some explicit cases I've tried of $A,B\in \mathbb{Z}$ appears to give me a flat $D\to S$.
For instance, if $E/\mathbb{Z}[6^{-1}]$ is the curve with  $(A,B)=(0,1)$ (which reduces to a supersingular curve at $p=5$), I find
$$h= 729z^{8}-1350x^4z^4+360x^6z^2+5x^8.$$
