Let $y^2-x^3 = c$ be Bachet's equation and pretend $(x,y)$ is a solution.
The tangent at $(x,y)$ of Bachet's curve is going to intersect it in a unique new point whose coordinates are supposed to yield the so-called "duplication formula":
$(\frac{x^4-8cx}{4y^2},\frac{-x^6-20cx^3+8c^2}{8y^3})$
I have not been able to re-derive this formula no matter how hard I try. It just won't work. Attempt:
Let $f(x,y)=y^2-x^3-c$, then the curve is implicitly parametrized by $f(x,y)=0$. The gradient of $f$ at $(x,y)$, which is $(-3x^2,2y)$, is orthogonal to the curve at $(x,y)$.
Therefore, the equation of the tangent is $-3x^2(X-x)+2y(Y-y) = 0$
(I'm looking at the only line orthogonal to the gradient at $(x,y)$ going through $(x,y)$)
So, if $(X,Y)$ is the coordinates of the point of intersection I'm looking for, it should be the unique solution (different from $(x,y)$), of the following conditions:
$-3x^2(X-x)+2y(Y-y) = 0$ and $Y^2 - X^3 - c = 0$
with the assumption that, of course, by hypothesis, $y^2 - x^3 = c$
Is my reasoning sound until now?
Mathematica does not seem to agree with me, but it could be that the full-simplify function isn't sophisticated enough (it does not reduce to Bachet's formula)
FullSimplify[Solve[{2*y*(Y - y) - 3*x^2*(X - x) == 0,
Y^2 - X^3 - c == 0}, {X, Y}], y^2 - x^3 == c]
I'm not expecting anybody to actually carry out the tedious calculations, though if somebody can manage to make Mathematica or some other software spit out the formula on its own I would be happy as I am trying to program an algorithm which automatically takes a curve in and finds the formula
Your reasoning is correct. Substitute $Y = \dfrac{3x^2}{2y} X - \dfrac{x^2 - 2c}{2y}$ into $Y^2 = X^3 + c$ to obtain a cubic for $X$. Just don't brute force it. The key point is to observe that two of its roots are known, and equal to $x$ (do you see why?). Use a sum of roots theorem to get the third one.
BTW, to use the sum of roots, yo don't need the entire cubic, but just the coefficient at $X^2$, which is $\dfrac{9x^4}{4y^2}$. Now the root you are interested in is indeed $$\dfrac{9x^4}{4y^2} - 2x = \dfrac {9x^4 - 8xy^2}{4y^2} = \dfrac{9x^4 - 8x(x^3 + c)}{4y^2} = \dfrac{x^4 - 8cx}{4y^2}$$
Proceed for $y$.