$\DeclareMathOperator{Pic}{Pic}$The dual isogeny can be shown to exist abstractly, but it can also be defined explicitly using the Picard group of classes of degree zero divisors, see Silverman's Arithmetic of Elliptic Curves III.6. This goes as follows:
If
$ \phi \colon E_1 \to E_2$
is an isogeny then
$$\hat\phi \colon E_2 \to \Pic^0(E_2) \to \Pic^0(E_1) \to E_1$$
is the dual isogeny.
What are all the maps here?
Well, $E \to \Pic^0(E)$ is simply $P \mapsto (P) - (0)$ which has the inverse map $(P) \mapsto P.$ We use the group structure on both sides to find the images of other divisors.
The map $\Pic^0(E_2) \to \Pic^0(E_1)$ is the pullback $\phi^*$ which sends $(P) \mapsto \sum_{Q \in \phi^{-1}(P)} (Q).$
So this is a pretty explicit recipe: to find what happens to $P\in E_2$ we just have to find the preimages of $P$ and of $0$ under $\phi,$ and sum them all in $E_1.$
Let's do an example. If I take a look at the LMFDB, I find some not too bad looking curves with a 3-isogeny between them. The curves are
$$E_1 = 14.a5\colon y^2 + x y + y = x^{3} - x$$
and
$$E_2 = 14.a6\colon y^2 + x y + y = x^{3} + 4 x - 6.$$
The 3-isogeny $\phi \colon E_1 \to E_2$ we'll use is given by
$$
\left(\frac{x^3 - x + 1}{x^2},\frac{x^3y + x^2 + xy - x - 2y - 1}{x^3}\right).
$$
(I used SageMath to get the equation for this.)
Let's find what $\hat\phi$ of the point $(1,-1) \in E_2(\mathbf{Q})$ is.
First we need to find $\phi^{-1}(1,-1),$ so we have to solve
$$
\frac{x^3 - x + 1}{x^2} = 1,
$$
$$
\frac{x^3y + x^2 + xy - x - 2y - 1}{x^3} = -1,
$$
which simplifies to
$$
x^3 - x^2 - x + 1 = 0,
$$$$
y(x^2 + 2x - 3) + (2x^2 - 2) = 0.
$$
The first polynomial factors as $(x-1)^2(x+1),$ so the $x$-coordinates of the preimage are $1$ and $-1.$ Labelling these roots as $\alpha_i,$ the preimage points are
$$\left(\alpha_i, -\frac{2\alpha_i^2 - 2}{\alpha_i^2 + 2\alpha_i -3}\right)=\left(\alpha_i, -\frac{2\alpha_i + 2}{\alpha_i+3}\right) \text{
for }\alpha_i \ne 1$$
When $x = 1$ the equation for the $y$-coordinates degenerates. To find the $y$-coordinates, we have to look at the original curve equation $$y^2 + 1\cdot y + y = 1^3 -1 =0.$$ So we have $y=0$ and $y=-2,$ giving $(1,0),(1,-2),$ and $(-1,0)$ as the preimages.
On the other hand, the non-zero preimages of the point $0 \in E_2(\mathbf{Q})$ are the points where $x = 0$ in order for the rational function defining $\phi$ to give $\infty.$ The y-coordinates of the two points on $E_1$ where $x = 0$ are the solutions of $y^2 + y = 0,$ i.e. $0$ and $-1$. We also have the point at infinity on $E_1.$
So $\phi^*((P) - (0)) = ((1, 0)) + ((1, -2)) + ((-1, 0)) - (0) - ((0,0)) - ((0,-1));$ hence, the point $\hat \phi (P)$ is just this sum evaluated on $E_1$, which is simply $(-1, 0)$ because points with the same $x$-coordinate cancel.
Most of the steps here were just fiddling around with functions. Even the addition at the end is just about complicated rational functions, but it's still surprising that one can in fact find nice enough rational functions for the dual. In this case, SageMath tells me that
$$
\hat \phi = \left(\frac{(x - 2)(x^{2} + 2 x + 13)}{(3 x + 1)^2}, \frac{(x + 5)(- x^{3} + x^{2} y - 8 x^{2} - 4 x y + 9 x + 11 y + 8)}{(3 x + 1)^3}\right).
$$
One final point is that the fundamental property of the dual isogeny allows us to do away with much of this. If $\hat\phi\circ\phi = [\deg \phi]$ then to find $\hat \phi (P)$ we can find $Q$ with $\phi(Q) = P$ so that $$\hat\phi(P) = \hat\phi(\phi(Q)) = [\deg \phi] Q.$$
In the above example, this means we could have stopped once we found the point $(-1,0)\in E_1(\mathbf Q)$ and tripled it to find the answer instead as we have a degree 3 isogeny. Since letting $x= -1$ in the equation for $E_1$ reduces to $y^2 = 0,$ the line $ x= -1$ is tangent to $E_1$ and $(-1,0)$ is 2-torsion; so, $$3\cdot(-1,0) = (-1,0)$$ just as before.
EllipticCurveIsogeny
class, and refers the Reader to Sage documentation for more details. – hardmath Dec 30 '17 at 12:57