A standard method to find such an equation (in principle) is the following, shown to me by Tate, but maybe due to Mordell(?). If $P_0$ is a point of order at least $4$ on $E$ (infinite order is allowed), then after a change of coordinates, we can find an equation for $E$ of the form
$$ E : y^2 + u x y + v y = x^3 + v x^2 \quad\text{with}\quad P_0=(0,0). $$
(It's a nice exercise to check this assertion.) So now we start with this equation and point, and we want $P_0$ to have order $N$. If $N$ is odd, say $N=2n+1$, then use the group law on $E$ to compute $x(nP_0)$ and $x\bigl((n+1)P_0\bigr)$. These will each be rational functions of $u$ and $v$; i.e., they're in $\mathbb Q(u,v)$. Setting
$$ x(nP_0)=x\bigl((n+1)P_0\bigr)$$
and clearing the denominators will give you a polynomial in $(u,v)$ whose vanishing is an affine plane model for $X_1(N)$ if $N$ is prime. However it tends to be rather singular at points where the discriminant of $E$ vanishes. And if $N$ is composite, then the polynomial will factor, i.e., the curve will have components for each divisor of $N$ greater than $3$. Finally, if $N=2n$ is even, you can do the same thing with the decomposition $N=(n+1)+(n-1)$, or you could set the numerator of $(2y+ux+v)(nP_0)$ equal to $0$, which is another way of forcing $2nP_0$ to be $O$.