Solution 1:

As you observe, the homogeneous coordinate $Y$ can never be zero on a $2$-torsion point, so we may work in the affine patch with coordinates $x = X/Y$, $z = Z/Y$, where the elliptic curve has equation $z^2 + z = x^3$, and origin (for the group law) equal to the point $(0,0)$.

In these coordinates, the map $[-1]$ is given by $(x,z) \mapsto (-\dfrac{x}{1+z},-\dfrac{z}{1+z})$. (This is not defined when $z = -1$, because $[-1]$ takes $(0,-1)$ to a point which lies at infinity in these coordinates, so we can remove the line $z = -1$ if you like.)

Thus, in these coordinates, $E[2]$ is cut out by the equations $z+z^2 - x^3 = 2x+xz = 2z + z^2 = 0$, just as you indicate, which (as you compute) reduce to the equations $x^4 - 2x = 0$ and $z = -x^3$.

But now you can just use the usual chord tangent law for addition:

Given points $P_1 = (x_1,z_1)$ and $P_2 = (x_2,z_2)$, let $m$ be the gradient of the line joining $P_1$ and $P_2$. The third point of intersection of this line with $E$ then has $x$-coordinate equal to $m^2 - x_1 - x_2$.

The gradient $m$ is given by the formula

$$\dfrac{z_1-z_2}{x_1 - x_2} = \dfrac{(z_1-z_2)(x_1^2 + x_1 x_2 + x_2^2)}{x_1^3 - x_2^3} = \dfrac{(z_1-z_2)(x_1^2 + x_1 x_2 + x_2^2)}{z_1 - z_2 + z_1^2 - z_2^2} = \dfrac{x_1^2 + x_1 x_2 + x_2^2}{1 + z_1 + z_2}.$$

Now we can restrict attention to $E[2]$, where $z = -x^3$, to obtain the addition formula (written in terms of $x$-coordinates alone, because points in $E[2]$ are determined by their $x$-coordinates)

$$(x_1,x_2) \mapsto \dfrac{(x_1^2 + x_1 x_2 + x_2^2)^2}{(1 -x_1^3 - x_2^3)^2} - x_1 - x_2.$$

Of course you can simplify this using the relations $x_1^4 = 2 x_1$ and $x_2^4 = 2 x_2$; I will leave that to you.

Added: In light of a comment below, maybe I should point out that this formula does answer the OP's question!

Namely, $x_1$ and $x_2$ are just more convenient names for the elements $x\otimes 1$ and $1 \otimes x$ of the tensor product. The above formula can then be rephrased by saying that the comultiplication is given by

$$x \mapsto \dfrac{(x_1^2 + x_1 x_2 + x_2^2)^2}{(1 - x_1^3 - x_2^3)^2} - x_1 - x_2.$$

This formula makes sense, because $(1- x_1^3 - x_2^3)$ is invertible in the ring

$$\mathbb Z_2[x_1,x_2]/(x_1^4 - 2 x_1, x_2^4 - 2 x_2) = \mathbb Z_2 [x]/(x^4 - 2 x) \otimes \mathbb Z_2[x]/(x^4 - 2x).$$

One more thing: the elliptic curve $E$ has good reduction over $\mathbb Z[1/3]$, and the above formula actually computes $E[2]$ as a finite flat group scheme over $\mathbb Z[1/3]$. (I.e. we can replace $\mathbb Z_2$ by $\mathbb Z[1/3]$ in the discussion.)