Homogenization refers to the process of going from the "affine" to the "projective" plane, by adding "points at infinity". It's easier to work over the projective plane because there aren't any "missing points"; for example, in the projective plane, every 2 distinct lines intersect at one point (instead of having to make an exception for parallel lines, as in the affine case). In effect, for each class of parallel lines, we add a point at infinity to the plane, and 2 parallel lines will intersect at the point at infinity corresponding to their class; the set of all points at infinity from a line at infinity.
For the study of elliptic curves, adding these missing points is necessary because otherwise, the points of the curve will not form a group (you'll have some missing points, like the point $O$ in the article you linked to).
Here's how the process works. I'll just show how it works in the plane; the same thing works in any number of dimensions. In the usual ("affine") plane, a point is given by coordinates $(x,y)$. In the projective planes, the coordinates will be written $(x:y:z)$. However, unlike the affine case, not all sets of coordinates will be distinct points; in the projective plane,
- (0:0:0) is not a point.
- If 2 coordinates differ by a scalar multiple, we consider them the same point: $(x:y:z) \sim (ax:ay:az)$. Here $a \ne 0$, and $x,y,z$ cannot all be 0, because of condition 1.
Just like in the affine case, we can write equations of curves using the projective coordinates $(x:y:z)$. However, we have to be careful, because of the equivalence condition imposed by condition 2 above. We cannot write a projective equation like $x+y+z=1$, because if we scale everything by a scalar $a$, we would get $ax+ay+az=1$, which is not equivalent to $x+y+z=1$. Indeed, condition 2 limits us to equations of the form $f(x,y,z)=0$ where $f$ is a homogeneous polynomial, meaning all terms of the polynomial have the same degree. For example, $y^2z - x^3 - 2xz^2 - z^3 = 0$ is a legitimate equation for a curve in the projective plane, because if we scale each coordinate by a scalar $a$, all the terms are multiplied by $a^3$, which can then be cancelled out (because $a \ne 0$ and the right-hand side is 0!). You can "homogenize" an affine equation in coordinates $x,y$ by multiplying each term by an appropriate power of $z$. That's what I did to get the equation above; starting with the affine elliptic curve equation $y^2 - x^3 - 2x - 1 = 0$, I homogenized each term to get the projective equation.
Because of condition 2, whenever we have a point $(x:y:z)$ and $z \ne 0$, we can normalize the $z$ coordinate to 1: $(x:y:z) \sim (x/z, y/z, 1)$. So among the points with non-zero $z$ coordinate, we effectively have only 2 coordinates to worry about, and we are back in the affine case. So the affine plane is contained in the projective plane: $(x,y) \mapsto (x:y:1)$. But the projective plane also contains an additional set of points with $z$-coordinate 0. These are the points at infinity.
In the affine plane, the lines $x+y=1$ and $x+y=2$ are parallel. Homogenizing, we get $x+y=z$ and $x+y=2z$ in the projective plane, where they intersect at a point at infinity (1:-1:0). (You can check that any other solution to this pair of equations is equivalent, via condition 2, to (1:-1:0), so there is only one intersection point in the projective plane.) The point (1:-1:0) is the point at infinity corresponding to (affine) lines of slope 1. We added one point an infinity for each possible slope (including $\infty$ for vertical lines); (1:m:0) for lines of slope $m$, and (0:1:0) for vertical lines.
In the elliptic curve case, homogenizing adds the point (0:1:0) to the curve, which is the point $O$ in the article linked in your original post; this point is usually taken to be the identity of the group. As noted above, (0:1:0) is the point at infinity corresponding to vertical lines. If you take an affine elliptic curve such as $y^2 - x^3 - 2x - 1 = 0$ and try to add the 2 points $(x,y)$ and $(x,-y)$, when you draw the line through them, they don't intersect the curve in a third point. But in the projective plane, they do; $(x,y)$ and $(x,-y)$ lie on a vertical line and therefore the line through them contains the point (0:1:0).
The situation with rational functions is similar. I'll switch down to 1 dimension since that illustrates the point more clearly. In the affine line, with a single coordinate $x$, I can define the rational function $1/x$. This has a pole at $x=0$ but no zeroes. However, on the projective line (with projective coordinates $(x:y)$), I homogenize to the function $y/x$. Now there is a pole (0:1) as before but there is also the zero (1:0). There are the same number of zeroes as poles. This is true for every rational function (not identically zero) on the projective line. It is also true for elliptic curves and indeed for every projective curve.
Note: This answer may give the impression that the $z$-coordinate is somehow "special" in the projective plane. This isn't really the case; it's an artifact coming from the way I chose to embed the affine plane into the projective plane, via $(x,y) \mapsto (x:y:1)$. This is convenient for relating the projective plane back to the familiar case of the affine plane, but it's just a choice, sort of like choosing a basis for a vector space; I could also choose either of the other coordinates to be the special one, or even something crazy like $(x,y) \mapsto (x+y+1:2x+3y+4:4x+5y+7)$ (not that one would actually want do something like that, which would be like deliberately choosing a weird basis for a vector space).
Instead of concentrating on calculating the divisor of $y$, calculate first the divisor of $x-e_i$, for $i=1,2,3$.
First homogenize to $Y^2Z=(X-e_1Z)(X-e_2Z)(X-e_3Z)$. Let us calculate the divisor of $x-e_1$ (the divisors of $x-e_2$ and $x-e_3$ are calculated in the same manner). The original function $x-e_1$ is the function $(X-e_1Z)/Z$ in projective coordinates. But, from the equation of the curve we see that
$$\frac{X-e_1Z}{Z} = \frac{Y^2}{(X-e_2Z)(X-e_3Z)}.$$
Thus, it is clear now that the function $(X-e_1Z)/Z$ has a double pole at $[X,Y,Z]=[0,1,0]$ and a double zero at $[e_1,0,1]$ because $e_1\neq e_2$, and $e_1\neq e_3$, as the curve is non-singular. Thus,
$$\operatorname{div}((X-e_iZ)/Z) = \operatorname{div}(x-e_i) = 2P_i - 2\infty, $$
where $P_i = [e_i,0,1]$ and $\infty=[0,1,0]$. Hence,
$$ \operatorname{div}(y^2)=\operatorname{div}(x-e_1) +\operatorname{div}(x-e_2) +\operatorname{div}(x-e_3)=2P_1+2P_2+2P_3-6\infty$$
and
$$ \operatorname{div}(y)=P_1+P_2+P_3-3\infty.$$
Edit to add: Let me work this out in more detail. Let $K$ be a field, let $\overline{K}$ be a fixed algebraic closure, let $E/K$ be an elliptic curve given by
$$ZY^2=X^3+AXZ^2+BZ^3=(X-e_1Z)(X-e_2Z)(X-e_3Z),$$ with $e_i\in K$ being distinct (notice the short Weierstrass equation implies $e_1+e_2+e_3=0$). Let $\overline{K}[E]$ be the function field of $E$. For $P\in E(\overline{K})$, we denote the ideal of functions in $\overline{K}[E]$ that vanish on $P$ by $M_P$ (strictly speaking we should define $M_P$ to be the maximal ideal in the localization $\overline{K}[E]_P$). For $f\in M_P$, we define
$$\operatorname{ord}_P(f)=\sup\{d\geq 1: f\in M_P^d\}.$$
Some facts:
Let $P=P_1=[e_1,0,1]$. Then, $X-e_1Z$ and $Y$ belong to $M_{P}\subseteq \overline{K}[E]$, and moreover, they generate $M_{P}$ (because $\langle X-e_1Z,Y\rangle$ is maximal in $\overline{K}[E]$). But, in fact,
$$(e_1-e_2)(e_1-e_3)(X-e_1Z)Z^2=ZY^2-(X-e_1Z)^3-(2e_1-e_2-e_3)(X-e_1Z)^2Z.$$
Since $Z$ does not vanish at $P=P_1=[e_1,0,1]$, and $(e_1-e_2)(e_1-e_3)$ is a non-zero constant, it follows that $X-e_1Z$ belongs to $M_P^2$. So $M_P/M_P^2$ is generated solely by $Y$ and, in particular, $\operatorname{ord}_P(Y)=1$ (note that if $Y$ was also in $M_P^2$, then $M_P/M_P^2$ would vanish, but this is impossible because $E$ is smooth and $\dim_{\overline{K}}(M_P/M_P^2)=\dim E=1$). The equation above also says that
$$\operatorname{ord}_P(X-e_1Z)=\min\{2\operatorname{ord}_P(Y),2\operatorname{ord}_P(X-e_1Z),3\operatorname{ord}_P(X-e_1Z)\},$$
and since $\operatorname{ord}_P(X-e_1Z)>0$, the only non-contradictory statement is that $\operatorname{ord}_P(X-e_1Z)=2\operatorname{ord}_P(Y)=2$.
Similarly, one deduces that at $P_i$, for any $i=1,2$, or $3$, the function $X-e_iZ$ vanishes to order $2$, and $Y$ vanishes to order $1$.
Now let $P=\infty=[0,1,0]$. The ideal $M_\infty$ is generated by $X$ and $Z$. It follows that $X-e_iZ\in M_\infty$ for any $i=1,2,3$. But
$$ZY^2 = (X-e_1Z)(X-e_2Z)(X-e_3Z)$$
and the fact that $Y\not\in M_\infty$, implies that $Z$ itself belongs to $M_\infty^3$. We conclude that $M_\infty/M_\infty^2$ is generated by $X$, and so $\operatorname{ord}_\infty(X)=1$, and this implies that $\operatorname{ord}_\infty(Z)=3$, and $$\operatorname{ord}_\infty(X-e_iZ)=\min\{\operatorname{ord}_\infty(X),\operatorname{ord}_\infty(Z)\}=\min\{1,3\}=1,$$
for $i=1,2,3$.
Hence, $\operatorname{ord}_\infty((X-e_iZ)/Z)=1-3=-2$, so $(X-e_iZ)/Z$ has a pole of order $2$ at $\infty$. Also, $\operatorname{ord}_\infty(Y/Z)=0-3=-3$, i.e., $Y/Z$ has a pole of order $3$ at $\infty$.
Putting it all together in terms of divisors we find:
$$\operatorname{div}(Y/Z)=(P_1)+(P_2)+(P_3)-3(\infty), \text{ and } \operatorname{div}((X-e_iZ)/Z) = 2(P_i)-2(\infty).$$
Best Answer
This fairly old-fashioned argument is the way I look at the situation. Homogenize $y^2=x^3+ax+b$ to $Y^2Z=X^3+aXZ^2+Z^3$, then dehomogenize by setting $Y=1$ to get $\zeta = \xi^3 +a\xi\zeta^2+\zeta^3$. The point you’re interested in is now $(0,0)$.
What does the curve look like there? Clearly $\xi$ is a local uniformizer, and $\zeta$ has a triple zero there. If you don’t see that right away, use the equation relating the two letters to see that $\zeta$ expands as a power series that starts $\zeta = \xi^3 + a\xi^7 +\cdots$. Now, what are the original $x$ and $y$ in terms of $\xi$ and $\zeta$? Yes: $x=\xi/\zeta$ and $y=1/\zeta$, and there you are.