Deriving Bachet's duplication formula

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$$.