Why does the diophantine equation $x^2+x+1=7^y$ have no integer solutions? This following Problem is from Pell equation chapters exercise

Let $y>3$ positive integer numbers, show that following diophantine equation
  $$x^2+x+1=7^y\tag{1}$$ has no integer solutions.

I tried write the equation
$$(2x+1)^2+3=4\cdot 7^y$$
if $y=2k$ then we have
$$(2\cdot 7^k+2x+1)(2\cdot 7^k-2x-1)=3$$
this case has no integer.
But $y$ is odd number, How to prove equation  (1)  has no  integer solutions for $x,y (y>3)$? Any help would be appreciated.
 A: Sketch of My Idea:
Let $\omega:=\frac{-1+\sqrt{-3}}{2}$ and $\bar{\omega}:=\frac{-1-\sqrt{-3}}{2}$.  Then, $7=(2-\omega)(2-\bar{\omega})$.  Hence, 
$$(x-\omega)(x-\bar{\omega})=x^2+x+1=7^y=(2-\omega)^y(2-\bar{\omega})^y\,.$$
Since $7\nmid x-\omega$, we may assume without loss of generality that $x-\omega=u(2-\omega)^y$, where $u$ is a unit (i.e., $\pm1$, $\pm\omega$, and $\pm\bar{\omega}$).  Then, $x-\bar{\omega}=\bar{u}(2-\bar{\omega})^y$, where $\bar{u}$ is the complex conjugate of $u$.  Thus,
$$-1=\frac{u(2-\omega)^y-\bar{u}(2-\bar{\omega})^y}{\omega-\bar{\omega}}\,.$$
Now, define $a_n(u)$ to be $\frac{u(2-\omega)^n-\bar{u}(2-\bar{\omega})^n}{\omega-\bar{\omega}}$ for $n\in\mathbb{Z}$.  You have that $a_n(u)=5a_{n-1}(u)-7a_{n-2}(u)$ for all $n\in\mathbb{Z}$.  Then, investigate each $u$.  You should find that, only for $n=0$, $n=1$, or $n=3$, there exists a unit $u$ such that $a_{n}(u)=-1$.  Now, since $a_n(-u)=-a_n(u)$, you may only look for solutions to $a_n(u)=\pm 1$ with $u \in \{1,\omega,\bar{\omega}\}$.
A: A sketch of my thoughts:
Let $\nu_7(n)=\max\{m\in\mathbb{N}:7^m\mid n\}$. A good idea may be to prove that if $\nu_7(x^2+x+1)=\nu_7(x^3-1)-\nu_7(x-1)=y>3$ then $x$ has to be large, say $x\geq 7^{y-1}$. In such a case, however, $x^2+x+1$ is too big to be just $7^y$ and it must have some other prime factor.
By this way, the problem boils down to finding (or, at least, lower-bounding) the two elements of order three in $G=\left(\mathbb{Z}_{/(7^y\mathbb{Z})}\right)^*$, that is a cyclic group with $o(G)=6\cdot 7^{y-1}$.
If $y=1$, that elements are $2$ and $4$. If $y=2$, that elements are $18$ and $30$. If $y=3$, that elements are $18$ and $18^2$. In general, we may compute such elements by solving $z^2+3\equiv 0\pmod{7^y}$. If $y=4$, such elements are $1047$ and $1353$. If $y=5$, such elements are $1353$ and $15453$.
We just need to find a pattern, or an alternative reason for which $4\cdot 7^y-3$ cannot be a square for $y>3$. Obviously if $y$ is even $4\cdot 7^y-3$ is too close to a square to be a square itself, so we may assume that $y$ is odd.
A: While not exactly solving your question, we can show that it is possible to find all solutions to a more general class of equations which includes the one asked for here. This shows that the equation asked for has only finitely many solutions and gives a method of finding them all - so, it shows that the question can be solved by performing some (rather tedious) calculations.
Letting $Q\in\mathbb Q[X]$ be a rational quadratic polynomial and $\mathcal P$ be a finite set of primes, then there are only finitely many solutions to
\begin{align}
Q(x)\in\mathcal P^*&&{\rm(1)}
\end{align}
for $x\in\mathbb Z$, where $\mathcal P^*$ is the set of integers whose prime factors are all in $\mathcal P$. Furthermore, the set of solutions is effectively computable.
In the question asked here, we have $Q=X^2+X+1$ and $\mathcal P=\{7\}$, and want to show that the only solutions are $x=-19,-3,-1,0,2,18$, for which $Q(x)\in\{1,7,7^3\}$. Now, this answer does not quite answer your question -- as I am going to apply an effective theorem on diophantine approximation to quadratic irrationals, although I haven't yet worked out how to compute the solutions. The theorem says that an algorithm exists, although a closer reading of the proof should give the algorithm (whether it is simple enough to be easily applied is another matter though).
To prove that (1) has finitely many solutions, we can complete the square to write
$$
a^{-1}\left(L(x)^2-b\right)\in\mathcal P^*
$$
for a degree one polynomial $L\in\mathbb Z[X]$ and fixed integers $a,b$. This can then be written as
$$
L(x)^2-b=au^2v
$$
for nonnegative $u\in\mathcal P^*$ and $v$ one of the finitely many squarefree elements of $\mathcal P^*$. In the example given, we have the case $v=1$ which is solved in the question, and $v=7$ which is asked for. We can rearrange as
$$
\left(L(x)+u\sqrt{av}\right)\left(L(x)-u\sqrt{av}\right)=b.
$$
From this, we obtain the inequality
$$
0 < \lvert u\sqrt{av}-w\rvert u\sqrt{av}\le\lvert b\rvert,
$$
where I have set $w=\lvert L(x)\rvert$. Rearranging, we express this as a rational approximation to the algebraic number $\sqrt{av}$,
\begin{align}
0 < \left\lvert\sqrt{av}-\frac wu\right\rvert\le \frac{\lvert b\rvert}{\sqrt{av}}u^{-2}.&&{\rm(2)}
\end{align}
Now, a result on diophantine approximation says that if $\xi$ is a real quadratic number then there exists effectively computable positive reals $C,\kappa > 0$ such that
$$
\left\lvert\xi-\frac wu\right\rvert\ge Cu^{-2+\kappa}
$$
for all integers $u,w$ with $u\in\mathcal P^*$. See the statement and proof of Theorem 1.2 in Effective Results For Restricted Rational Approximation To Quadratic Irrationals by Bennet and Bugeaud (link). Taking $\xi=\sqrt{av}$ and comparing with (2),
$$
Cu^{-2+\kappa}\le\frac{\lvert b\rvert}{\sqrt{av}}u^{-2}.
$$
Hence,
$$
u\le\left(\frac{\lvert b\rvert}{C\sqrt{av}}\right)^{1/\kappa}.
$$
So, we have a computable upper bound for $u$, and there are only finitely many solutions to (1).
A: An elementary solution [albeit not using Pellian methods] is reasonably simple to obtain.
As you show, there is no solution when $y$ is even. Now assume $y$ is odd. Observe that $7 \mid (x^2+x+1)$ implies $x \equiv 2\!\pmod{7}$ or $x \equiv 4\!\pmod{7}$. 
In the first case, we substitute $x=7z+2$ into the equation and simplify to obtain
$$z(7z+5)=\bigl(7^\frac{y-1}{2}-1\bigr)\bigl(7^\frac{y-1}{2}+1\bigr).$$
Now our assumption of a positive solution with $y \ge 3$ implies there exist positive integers $a,b,c,d$ such that
\begin{align}
z &= ab, & 7^\frac{y-1}{2}-1 &= ac,  \\
7z+5 &= cd, & 7^\frac{y-1}{2}+1 &= bd.
\end{align}
The conclusion follows fairly directly via algebraic manipulations.
The proof in the second case is nearly identical.

EDIT:

Here's a hint to get most of the way through the proof… Writing $k=(y-1)/2$, we have
$$
\frac{z}{7^k+1} = \frac{a}{d} = \frac{7^k-1}{7z+5}.
$$
Solving each relation for $7^k$ gives
$$
7^k = \frac{dz-a}{a} = \frac{7az+5a+d}{d},
$$
and solving the second relation for $z$ yields
$$
z = \frac{a(2d+5a)}{d^2-7a^2}.
$$
Since $z$ is an integer by hypothesis, there must be a prime $p$ dividing $\gcd(d,a)$, or $(d^2-7a^2) \mid (2d+5a)$, from which we can quickly deduce $p \mid 3a$.
A: We work in the Eisenstein integers. The equation then becomes$$(x + 1 + \omega)(x - \omega) = 6^y.$$Suppose $d\, |\, (x + 1 + \omega)$ and $\text{ }d\,|\,(x - \omega)$.Then we have $d\,|\,(1 + 2\omega)$, so $N(d)\,|\,3$. However, we also have $N(d)\,|\,7^{2y}$, so $N(d) = 1$. On the other hand, $7 = (-2 - 3\omega)(1 + 3\omega)$ and $-2 - 3\omega$, $1 + 3\omega$ are relatively prime. Thus, it follows that $x + 1 + \omega$ and $x - \omega$ are $e(-2-3\omega)^y$ and $f(1 + 3\omega)^y$ in some order, where $e$ and $f$ are units $\pm1$, $\pm\omega$, $\pm1\pm\omega$ with $ef = 1$. We just need to show that if$$a + b\omega = f(1 + 3\omega)^y,$$then $|b|$ can never be $1$, which would solve the problem. Indeed, let$$A = 1 + 3^3\binom{y}{3} + \dots,\text{ }B = 3\binom{y}{1} + \dots,\text{ }C = 3^2\binom{y}{2} + \dots.$$We have$$a + b\omega = f(A - C + \omega(B - C)),$$so we want to show none of $|A - B|$, $|B- C|$, $|C - A|$ are one. However, this is a relatively easy task for $y > 3$, so we are done.
