Let $k$ be a natural number . Then $3k+1$ , $4k+1$ and $6k+1$ cannot all be square numbers.

Complete answer from the literature:

I found a list of solved Pell systems in Szalay, Appendix 4, page 84.

Kiran Kedlaya published Solving Constrained Pell Equations in Mathematics of Computation, Volume 67, April 1998, pages 833-842.

As an example of his method, on page 840 he does the harder part of the Lucas problem (bottom page 238), namely

Solving $n+1 = 2 x^2,$ $2n+1 = 3 y^2,$ $n=z^2$

with the conclusion

Possible values of $n:$ $\{1\}$

This fits the problem above by taking $x^2 = 3k+1,$ $y^2 = 4k+1,$ $z^2 = 6k+1$ which gives the system $$ 2 x^2 - z^2 = 1,$$ $$ 3 y^2 - 2 z^2 = 1, $$ with the conclusion that we can only have $z^2 = 1,$ so $k=0.$

Kedlaya refers to his own solution as elementary, references bottom of page 838 to top of page 839. He mentions other articles with elementary methods; I think this is a case where "elementary" is in the eye of the beholder.

$\def\Q{\mathbf{Q}}$ $\def\Qbar{\overline{\Q}}$ $\def\Z{\mathbf{Z}}$ $\def\F{\mathbf{F}}$

I want to try to explain why there won't be any slick elementary argument in this case. The usual competition techniques involve a tricky use of congruences. A more advanced trick is to use infinite descent. However, in the end, this problem boils down to finding all integral points on an elliptic curve with positive rank (the curve $Y^2 = X^3 - 36 X$). There are some pretty standard techniques to do this - but none of them are elementary. The reason why elementary arguments (such as congruences) are doomed to fail is related to the fact that there are many solutions to this problem with $k \in \Q$.

Any odd square is $1 \mod 8$. Hence $k$ is even, and so we need to find $k$ such that $6k+1$, $8k+1$, and $12k + 1$ are squares.

As others have observed, can write down equations for this problem as $a^2 = 6k+1, b^2 = 8k+1$, and $c^2 = 12k+1$. This leads to the sytem equations $$C: 3 b^2 - 4 a^2 = -d^2, \quad c^2 - 2 a^2 = -d^2,$$ where $d = 1$. Here I introduced the variable $d$ to turn the affine equation into a projective one. The intersection of two quadrics in $\mathbf{P}^3$ has genus one, and so this equation is an elliptic curve (for some choice of base point, note that $[1;1;1;1] \in C(\Q)$.) It is a general theorem (of Siegel) that, given an elliptic curve has only finitely many integral points. On the other hand, most elementary solutions to these sort of problems typically prove the stronger claim that $E(\mathbf{Q})$, the set of rational points, is also finite. Proofs of this kind of statement were first found by Fermat, and use the method of infinite descent. For example, Fermat proved that there are no integral solutions to $x^4 + y^4 = z^4$ by considering the equation $x^4 + y^4 = z^2$ in rational numbers. After scaling, one can let $y = 1$, and then $z^2 = x^4 + 1$ is an affine open inside an elliptic curve. Unfortunately, as we shall see, this is not one of those examples.

For those who like to think of elliptic curves in terms of Weierstrass equations, it's easy to write down an elliptic curve $E$ isogenous to $C$ of the form $y^2$ is a cubic. Namely, if $6k + 1$, $8k+1$, and $12k+1$ are all squares, then so is their product, hence we also get a solution in integers (and hence in rational numbers) to the equation

$$E: z^2 = (6x+1)(8x+1)(12x+1),$$

Note that by replacing $k$ with $2k$ (which makes no difference to rational solutions), it makes the coefficient of $x^3$ equal to a square ($576 = 24^2$). And now, letting $z = 24y$, we can also write as

$$E: y^2 = (x + 1/6)(x + 1/8)(x + 1/12),$$

and the point $\infty$ at infinity is a rational point. There is a degree $2$ isogeny $C \rightarrow E$ given by $$(a,b,c) \mapsto \left(\frac{a - 1}{6},\frac{abc}{24} \right).$$ There will also be a dual isogeny $E \rightarrow C$ of degree $2$, but it will be a little messy to write down.

Descent: The main elementary way to find rational points on an elliptic curve $E$ is to do a $2$-descent. In the case when all the $2$-torsion is rational, this is pretty concrete. Explicitly, there is an injective map: $$E(\Q) \rightarrow \Q^{\times}/\Q^{\times 2} \times \Q^{\times}/\Q^{\times 2} \times \Q^{\times}/\Q^{\times 2}$$ which sends a point $P = [x,y]$ to the class: $$(x + 1/6, x + 1/8,x + 1/12).$$ Note that this class lies inside the subgroup $(a,b,c) \subset (\Q^{\times}/\Q^{\times 2})^3$ with $abc = 1$, so we can really just take the first two entries in $(\Q^{\times}/\Q^{\times 2})^2$. The key fact about this map is that the kernel is precisely $2 E(\Q)$, so the induced map from $E(\Q)/2 E(\Q)$ is injective. A slight subtlety --- this doesn't quite make sense for $2$-torsion points, because one of the entries is zero. But you can work out the correct image by noting that it should land inside the diagonal. Hence $$\infty \mapsto (1,1,1),$$ $$e_1 = [-1/6,0] \mapsto (*, -1/24,-1/12) \equiv (2,-6,-3) \in (\Q^{\times}/\Q^{\times 2})^3,$$ $$e_2 = [-1/8,0] \mapsto (1/24,*,1/24) \equiv (6,-1,-6) \in (\Q^{\times}/\Q^{\times 2})^3,$$ $$e_1 + e_2 = [-1/12,0] \mapsto (1/12,1/24,*) \equiv (3,6,2) \in (\Q^{\times}/\Q^{\times 2})^3.$$

On the other hand, $E$ also has an obvious rational point $P = [0,1]$. We find that $$P \mapsto (6,2,3) \in (\Q^{\times}/\Q^{\times 2})^3,$$ $$P + e_1 \mapsto (3,-3,1) \in (\Q^{\times}/\Q^{\times 2})^3, $$ $$P + e_2 \mapsto (1,-2,-2) \in (\Q^{\times}/\Q^{\times 2})^3,$$ $$P + e_3 \mapsto (2,3,6) \in (\Q^{\times}/\Q^{\times 2})^3.$$

Let's now make the following observation, using the fact that $\mathbf{Q}$ has unique factorization. If $(x + 1/6)$, $(x + 1/8)$, and $(x + 1/12)$ have a square product, then they are all squares up to possible common factors. (Any power of $p > 3$ in the numerator will have to be a square as well by an elementary computation). However, any common factor must divide $$\Delta = (1/6 - 1/8)^2(1/6 - 1/12)^2 (1/8 - 1/12)^2 = 2^{-16} \cdot 3^{-6}.$$ This implies that the image of $E(\mathbf{Q})/2 E(\mathbf{Q})$ lands in the smaller group $$(S^{\times}/S^{\times 2 } )^2,$$ Where $S \subset \mathbf{Q}^{\times}$ is generated by $2$, $3$, and $-1$. This is a group of order $8 \times 8 = 64$. We have computed a subgroup of order $8$ that lies in the image. If we are lucky, we can show that none of the other $56$ points lie in the image for local reasons, and this would imply that $E(\mathbf{Q})/2 E(\mathbf{Q}) = (\mathbf{Z}/2 \mathbf{Z})^3$. We don't need to check all $56$ cases, because the image is a group. And most of them are easy. Namely, we write: $$x + \frac{1}{6} = a u^2, x + \frac{1}{8} = b v^2, x + \frac{1}{12} = a b w^2,$$ or $$a u^2 - b v^2 = \frac{1}{24}, \qquad a u^2 - a b w^2 = \frac{1}{12},$$ where $(a,b) \in (S^{\times}/S^{\times 2 } )^2$. We can immediately rule out many examples: If $b > 0$ and $a < 0$ the first equation has no solutions. If $a < 0$ and $b < 0$, the second equation has no solutions. Hence $a > 0$. This leaves us with $32$ possibilities. Considering the equation $2$-adically rules out half of the remaining possibilities, and $3$-adic considerations reduce us by half again to end up with $8$. Thus we have shown:

$$E(\mathbf{Q}) \simeq \mathbf{Z} \oplus (\mathbf{Z}/2 \mathbf{Z})^2,$$ and $P = [0,1]$ is a generator. (Computing the torsion subgroup is easy: note that $|E(\F_{11})| = 12$ and $E(\F_{5}) = 8$, so the torsion subgroup has order dividing $4$.)

So we have now found (in some sense) all rational points on $E$. However, the answer is a little unexpected, because we had hoped to show that the original equation had no solutions by showing that there were no rational solutions. But there are many rational solutions!

Let's think of solutions to $E(\mathbf{Q})$ that actually come from $C(\mathbf{Q})$. We want $6x+1$, $8x +1$, and $12x+1$ to all be squares. Hence, we want the image of the point in $E$ in $(\mathbf{Q}^{\times}/\mathbf{Q}^{\times 2}$ to be $$(1/6,1/8,1/12) = (6,2,3).$$ Clearly $P$ is such a point, and this corresponds to $x = 0$ or $k = 0$. Any other point which maps to this must be equal to $P$ in $E(\mathbf{Q})/2 E(\mathbf{Q})$, and thus $P \mod 2 E(\mathbf{Q})$. Hence we deduce:

Rational points on $C(\mathbf{Q})$ correspond to points on $E(\mathbf{Q})$ of the form $2Q + P$, or equivalently, $n P$ for $n$ odd.

For example, $$3P = \left[ \frac{-140}{2209},\frac{-28083}{103823} \right],$$ and we correspondingly see that $$3 \left(\frac{-280}{2209}\right) + 1 = \left(\frac{37}{47}\right)^2, 4 \left(\frac{-280}{2209}\right) + 1 = \left(\frac{33}{47}\right)^2, 6 \left(\frac{-280}{2209}\right) + 1 = \left(\frac{23}{47}\right)^2$$ Taking $5P$ gives $$3 \left(\frac{369572280}{537451489}\right) + 1 = \left(\frac{40573}{23183} \right)^2, 4 \left(\frac{369572280}{537451489}\right) + 1 = \left(\frac{44897}{23183}\right)^2, 6 \left(\frac{369572280}{537451489}\right) + 1 = \left(\frac{52487}{23183}\right)^2,$$ and so on.

But what does this say about our original question? A weaker question than asking for solutions with integral $x$ is to ask for solutions with $72x \in \mathbf{Z}$. Then we can make the transformation $$X = \frac{x}{144} - \frac{1}{8}, \quad Y = \frac{y}{1728},$$ to get the equation $$D: Y^2 = X^3 - 6^2 X.$$ This is isomorphic to $E$ over $\Q$, but it is a different integral model. This is a familiar curve; it is the congruent number curve for $d = 6$. There is an integral solution related to the fact that $(3,4,5)$ is a right angled triangle with area $6$. But there are other integral solutions as well, namely: $$[0,0], [\pm 6, 0], [-3, \pm 9], [-2, \pm 8], [12, \pm 36], [18, \pm 72], [294, \pm 5040].$$ This is, in fact, the complete list. If you like, you can confirm this in magma as follows:

$$\begin{aligned} & \ \texttt{> E := EllipticCurve([0,0,0,-36,0]);} \\ & \ \texttt{> IntegralPoints(E);} \\ & \ \texttt{ [ (-6 : 0 : 1), (-3 : -9 : 1), (-2 : 8 : 1), (0 : 0 : 1),} \\ & \ \texttt{ (6 : 0 : 1), (12 : -36 : 1), (18 : 72 : 1), (294 : -5040 : 1) ] } \end{aligned}$$

The observatio that $E$ has many rational points and lots of $S$-integral points for $S \in \{2,3\}$ suggests that there won't be a quick solution to the original problem. There are fairly standard methods to (certifiably) compute all the integral points of an elliptic curve $D$ given the Mordell-Weil group (which $\texttt{magma}$ uses), but they are not particularly elementary.

I think this will work, but I don't quite have time today to fill in all details.

If $3k+1=x^2,4k+1=y^2,6k+1=z^2$, then $2x^2+z^2=3y^2$. Letting $a=x/y,b=z/y$, this gives $2a^2+b^2=3$. We can parametrize the whole space of solutions like one usually does for Pythagoream triples: $(1,1)$ is an obvious point on the curve. Now let $t$ be any real and let $Y=t(X-1)+1$ be a line through it. It intersects the curve $2a^2+b^2=3$ at point $(1,1)$ and other point $(a,b)$ satisfying $2a^2+b^2=3$ and $b=t(a-1)+1$. Substituting, we get $a^2(2+t^2)+a(2t-2t^2)+(t^2-2t-2)=0$. If we interpret this as a quadratic in $a$, $1$ is one solution, so by Viete's formulas, the other is $a=(t^2-2t-2)/(t^2+2)$. Then $b=(-t^2-4t+2)/(t^2+2)$. Now we look at $x/z=a/b=(t^2-2t-2)/(-t^2-4t+2)$. Letting $t=n/m$ be rational, this gives $x/z=(n^2-2nm-2m^2)/(-n^2-4nm+2m^2)$.

Clearly $x,z$ are relatively prime, since $2x^2-z^2=1$. I claim $n^2-2nm-2m^2,-n^2-4nm+2m^2$ have $\gcd$ dividing $6$. This is because their sum is $6nm$, and if $p\mid n$ and $p\mid n^2-2nm-2m^2$, then $p\mid m$ and vice versa.

(to be continued)