Instead of following blindly Wikipedia's formulas, it is best to understand how to calculate $P+Q$, or $P+P$, given an elliptic curve. Let us assume for simplicity that the curve is given by $E:y^2=x^3+Ax+B$, and $P,Q\in E$.
In order to find $P+Q$, first find the equation of the line $L$ through $P$ and $Q$, find the third point $R$ of intersection of $L$ and $E$. Then $P+Q+R=\mathcal{O}$, so that $R=-(P+Q)$. In this case $-(x_0,y_0)=(x_0,-y_0)$, so we can find $P+Q$ as $-R$.
In order to find $P+P$, first find the tangent line $L$ to $E$ at $P$, and find the third point $R$ of intersection of $L$ and $E$. Then, $2P+R=\mathcal{O}$, and so $2P=-R$.
In your case, $E: y^2=x^3+1$ and $P=(0,1)$. The tangent line is $y=1$, and it turns out that this line has a triple point of intersection with $E$, thus $R=P$. In particular, $2P=-R=-P$, and so $3P=\mathcal{O}$.
Another example with the same curve $E$: let $P=(2,3)$ and $Q=(0,1)$. The line through $P$ and $Q$ is $y=x+1$, and the third point of intersection is $R=(-1,0)$. Since $R=-R$, we obtain $P+Q=(-1,0)$.
Finally, here is a picture of $P=(2,3)$, $2P$, $3P$, $4P$, and $5P$, which shows that $6P=\mathcal{O}$.
When going from continuous fields ($\Bbb R$ and $\Bbb C$, and to a lesser extent $\Bbb Q$) to something like a finite field, you are right that the geometric intuition fails. However, nothing is stopping the algebra from working exactly as before (or, at least, almost exactly like before; as you pointed out, characteristic 2 and 3 mess with things when working with quadratic and cubic polynomials).
Say we have a finite field $K$. Then the line between $(a, b)$ and $(c, d)$ in $K^2$ is simply the set of $(x, y)\in K^2$ which satisfy a linear equation
$$
mx + ny + k = 0
$$
where $m, n, k$ are chosen so that the equation is satisfied for $(a, b)$ and $(c, d)$ (if the two points are distinct, then there is exactly one viable choice of $m, n, k$ up to scaling).
A tangent is a bit more tricky, but as we know from the continuous case, a line $mx + ny + k = 0$ which has a point $(a, b)$ in common with $E$ is tangent with $E$ at that point if the intersection at $(a, b)$ has degree greater than $1$. And the degree of the intersection can be found purely algebraically, for instance by looking at $K[x, y]/(mx+ny + k)\cong K[z]$ and considering whether $x^3 + Ax + B - y^2$ has a multiple root at $(a, b)$.
So the moral of the story is that even though the geometric interpretation doesn't make as much sense, we keep the geometric terms and define them algebraically instead, and much of the intuitive results from the continuous case carry right over to the discrete case.
Best Answer
You may be confusing the third point, say $S$, of intersection of $E$ and the line through $P$ and $Q$ with the point $P + Q = R$. Determining $S$ is just the first step in determining $R$. $R$ is, in turn, the third point of intersection of $E$ and line through $S$ and the zero point $O$. Wikipedia, as is traditional, chooses $O$ to be the point at infinity $\, 0:1:0 \,$ in the projective plane. With that choice for the zero point, the line through $S$ and $O$ is vertical, so, because $E$ is symmetric with respect to the $x$-axis (because $y$ appears only as $y^2$ in the equation for $E$), we just need to change the sign of the $y$-coordinate of $S$ to get $R$.
We now derive the Wikipedia formulas based on the two steps described above. The equation of the line through $P$ and $Q$ is $$y = y_p + \lambda (x - x_p).$$ Substitute that expression for $y$ into $y^2 = x^3 + ax + b$ to get $$x^3 - \lambda^2x^2 +(a + 2\lambda^2 x_p - 2\lambda y_p)x + b - (\lambda x_p - y_p)^2 = 0.$$ The three solutions to that cubic equation give the $x$-coordinates $x_p, x_q, x_s$ of the three points of intersection of the line with $E$. From Vieta's first formula, we see that the sum of those $x$-coordinates is $\lambda^2$ so that $x_p + x_q + x_s = \lambda^2$. When we reflect $S$ over the $x$-axis, the $x$-coordinate does not change, so $x_r = x_s$. Thus, $$x_r = \lambda^2 - x_p - x_q.$$ Using the equation of the line, \begin{align} y_s & = y_p + \lambda (x_s - x_p)\\ & = y_p + \lambda (x_r - x_p) \end{align} When we reflect $S$ over the $x$-axis, the sign of the $y$-coordinate changes, i.e., $y_r = -y_s$. Thus, $$y_r = \lambda (x_p - x_r) - y_p.$$