How to solve $(2k+1)(2k^2+1) = y^2$ over the positive integers

Question

How could I find a solution to $(2k+1)(2k^2+1) = y^2$ over positive integers?

WolframAlpha returns some solutions, but I was wondering if there was a way to find all of them, or prove that there are infinitely many. Here is the link to the Wolfram Alpha answer:

https://www.wolframalpha.com/input/?i=solve+%282k%2B1%29%282k%5E2%2B1%29+%3D+y%5E2+over+the+integers

This is a personal question inspired by the equation $(k+2)(4k+1) = y^2$ which is doable using elementary methods.

Answer 1

Here is another computer solution by using SAGEMATH (though we need to put it in standard form). We are indeed looking for integral points on an elliptic curve. First we transform our Diophantine equation $E_1:(2k+1)(2k^2+1) = y^2$ into the form $Y^2 = X^3 +a X^2 + b X +c$, where $a,b,c$ are integers (which SAGEMATH will compute integral points for us, assuming its algorithm is correct).

With the substitutions $k=X/64$ and $y=Y/256$, we obtain an elliptic curve in standard form with integral coefficients $E_2:Y^2 = X^3 + 32X^2+2048X+65536$. Note any integral solutions to our original problem $E_1$ will be an integral solution to this new equation $E_2$.

Using the following in SAGEMATH (which you can execute it here https://sagecell.sagemath.org/)

x,y = polygens(QQ,'x,y'); 
E = EllipticCurve(-y^2+x^3+32*x^2+2048*x+65536); 
E.integral_points(both_signs=True)

(One can also check this curve has rank 1 by using the command E.rank())

It shows the only integral solutions to $E_2$ are (in homogeneous coordinates $(X:Y:1)$)

[(-32 : 0 : 1),
 (-16 : -192 : 1),
 (-16 : 192 : 1),
 (0 : -256 : 1),
 (0 : 256 : 1),
 (64 : -768 : 1),
 (64 : 768 : 1),
 (160 : -2304 : 1),
 (160 : 2304 : 1),
 (768 : -21760 : 1),
 (768 : 21760 : 1)]

Note, by Siegel's Theorem, there can only be finitely many integral solutions (https://mathworld.wolfram.com/SiegelsTheorem.html).

Here we see there are only six solutions such that $X$ is a multiple of 64 and $Y$ is a multiple of 256, namely

 (0 : -256 : 1),
 (0 : 256 : 1),
 (64 : -768 : 1),
 (64 : 768 : 1),
 (768 : -21760 : 1),
 (768 : 21760 : 1)

Recovering for $k$ and $y$ we have the only six integral solutions offered in Wolfram Alpha: $(k,y) = (0,\pm 1), (1,\pm 3),(12,\pm 85)$.

I only have a superficial understanding of this, but I will direct you here for more information:

1. https://en.wikipedia.org/wiki/Elliptic_curve
2. Rational Points on Elliptic Curves, by Silverman.