To find the doubling map of the elliptic curve $y^2=x^3+1$ I want to verify that the doubling map of the elliptic curve $y^2=x^3+1$ is given by
\begin{align}
 P=(x,y) &\to 2P \\
 (x,y) &\mapsto \left(\frac{x^4-8x}{4x^3+4}, \frac{2x^6+40x^3}{8y^3} \right)
\end{align}

For calculation flexibility we denote $P=(x',y')$. We will find the tangent line at $P$ and make intersect with the curve to find $2P$.
Differentiating the curve with respect to $x$, we get $$\frac{dy}{dx}=\frac{3x^2}{2y}.$$
So the tangent line at $P$ is given by $$y-y'=\frac{3x'^2}{2y'}(x-x') \Rightarrow y=\frac{3x'^2x-3x'^3+2y'^2}{2y'}.$$
Substituting $y$ in the curve, we get
\begin{align}
&\left(\frac{3x'^2x-3x'^3+2y'^2}{2y'} \right)^2=x^3+1 \\
&\Rightarrow (3x'^2x-3x'^3+2y'^2)^2=4y'^2x^3+4y'^2 \\
&\Rightarrow 9x'^4x^2+9x'^6+4y'^4+18x'^5x+12x'^3y'^2-12x'^2y'^2x=4y'^2x^3+4y'^2 \\
& \Rightarrow 4y'^2x^3-9x'^4x^2+(12x'^2y'^2-18x'^5)x+(4y'^2-9x'^6-4y'^4-12x'^3y'^2)=0, \cdots (1)
\end{align}
The elliptic passes through $(0,1)$, and hence plugging it in $(1)$, we get
$$4y'^2-9x'^6-4y'^4-12x'^3y'^2=0, \cdots (2)$$
But now I am out of way, how to solve for $x', y'$ ?

what is the apropriate way to solve my question ?


I tried to check optically the claimed formula on an obvious example of choice, taking $Q=(0,1)$, since the two denominators are divisible by $x$, thus getting $2Q\overset?=(0,0)$. Not a point on the curve. What happens here?

My way of computing is as follows. Consider the point $P=(u,v)$ on the elliptic curve with affine equation $Y^2 =X^3 +1$. Then the tangent line to the given curve in $(u,v)$ is $2Y\; dY=3X^2\; dX$, which taken in $(u,v)$ is $2v(Y-v) = 3u^2(X-u)$. We isolate and replace $Y$ from this equation, it is $Y=mx+n$ with $$m=\frac{3u^2}{2v}\ ,$$
into $-Y^2 + X^3 +\dots$ to get $0 = X^3 -m^2X^2+\dots$, and we already know two roots (with multiplicity counted), which are $u,u$, so the third root, $x$, is given by (Vieta) $u+u+x=m^2$. So
$$
x = -2u+m^2 
= -2u + \frac {9u^4}{4v^2}
=  \frac {-8uv^2 + 9u^4}{4v^2}
=   \frac {-8u(u^3+1) + 9u^4}{4v^2}
=   \frac {u^4-8u}{4(u^3+1)}
\ .
$$
To get the $Y$-component $y'$ of $-2P=(x,y')$, we use
$$
\begin{aligned}
y'&=mx+n=m(-2u+m^2)+n=-2mu+m^3+(v-mu)
=v-3mu+m^3
\\
&=\frac 1{8v^3}(8v^4-36u^3v^2+27u^6)\text{ which is either}
\\
&=\frac1{8v^3}(8v^4-36(v^2-1)v^2+27(v^2-1)^2)
=\frac1{8v^3}(-v^4-18v^2+27)\text{ or}
\\
&=\frac1{8v^3}(8(u^3+1)^2-36u^3(u^3+1)+27u^6)
=\frac1{8v^3}(-u^6-20u^3\color{red}{+8})\ .
\end{aligned}
$$
This gives for $(x,y)=2P=-(x,y')=(x,y)$ the formula:
$$
2P=(x,y)=
\left(\ 
\frac {u^4-8u}{4(u^3+1)}\ ,\ 
\frac{u^6+20u^3\color{red}{-8}}{8v^3}\ \right)\ .
$$

Chech using the point $Q=(0,1)$: Applying the above formula we get $2Q=(0,-1)$, which is a point on the curve.

Sage check (using the implemented sage formula for the field $K$, which is the fraction field of the ring $R=\Bbb Q[U,V]/(V^2 -(U^3+1))$):
F.<U,V> = PolynomialRing(QQ)
R.<u,v> = F.quotient(V^2 - U^3 - 1)
K = FractionField(R)

E = EllipticCurve(K, [0, 1])
P = E.point([u, v])
print(f'P = {P.xy()}')
print(f'2P = {(2*P).xy()}')

This gives:
P = (u, v)
2P = ((u*v^2 - 9*u)/(4*v^2), (v^4 + 18*v^2 - 27)/(8*v^3))

