Describe in as much detail as possible the form of all positive integer solutions to the diophantine equation $2x^2+x = y+3y^2.$
Note that for any solution, $x-y, 3x+3y+1, 2x+2y+1$ must all be perfect squares. Also we cannot have $x\leq y$ as otherwise $2x^2+x\leq 2y^2+y < 3y^2+y$ so equality cannot hold. The equation is equivalent to $2x^2-3y^2=y-x$, or to $$3x^2-3y^2=y-x+x^2\Leftrightarrow (x-y)(3x+3y+1)=x^2 \quad(1).$$ If $x-y$ and $3x+3y+1$ are not coprime, then we can choose a prime p dividing both. So $p | x$ from the last equality. But then $p$ also divides $3x+3y+1+3(x-y) = 6x+1$, a contradiction. Also, ne has $$2x^2-3y^2 = y-x\Leftrightarrow (x-y)(2x+2y+1)=y^2\quad(2)$$ and so by similar reasoning to above, $x-y$ and $2x+2y+1$ are coprime (indeed if not, let $p$ divide both. Then $p|y$ and $p|(2x+2y+1-2(x-y)) = 4y+1,$ a contradiction). If $a$ and $b$ are two coprime numbers whose product is a square, then $a$ and $b$ must both be squares, since if one of them has an odd exponent of a prime in its prime factorization, then so will their product. Hence equations (1) and (2) show that $x-y,3x+3y+1,$ and $2x+2y+1$ must all be perfect squares.
Note that from the original equation, $x$ and $y$ are either both even, or $x$ is even and $y$ is odd (one can see this by considering separately the cases where $x$ and $y$ have the same or different parities and considering the parity of both sides of the equation). Write $$x-y=a^2,3x+3y+1=b^2, 2x+2y+1 = c^2\quad (0)$$ for positive integers $a,b,c$, where by (1), $ab = x$ and by (2), $ac = y$. Observe that $a < c < b.$ One can show that x and y are both composite using these two equations since $b > 1$ and if $a = 1$ so that $x=y+1$ and $b=x$, we have $3(y+1)+3y+1 = (y+1)^2 = y^2+2y+1\Rightarrow y^2-4y-3 = 0,$ which has no integer solutions for y. Also, $c>1$ so $y$ is composite. One can also deduce from the above equations that $c^2 = b^2- (x+y)=b^2-a(b+c),$ and so $c$ is a root of the quadratic equation $c^2+ac+b(a-b) = 0.$ We have by the quadratic formula that $c = \dfrac{-a\pm (2b-a)}{2} = b-a,$ because $c>0$. $a<c=b-a$ so $b>2a.$ $a>1$ as well. So as a specific example consider $a=2, b=5,c=3$. We have $x=y+4$. $a(b-c)=a^2, 3ab + 3ac+1 = b^2 =3a(b+c)+1 = 3a(2b-a)+1.$ So $b^2 - 6ab+3a^2-1=0\Rightarrow b = \dfrac{6a\pm \sqrt{36a^2 - 4(3a^2-1)}}{2}.$ $36a^2-4(3a^2-1) = 24a^2+4 = 4(6a^2+1)$ so $b=3a\pm \sqrt{6a^2+1}.$ Thus b is an integer if and only if $6a^2+1$ is a perfect square. Also $b>2a$ implies $3a \pm \sqrt{6a^2+1} > 2a.$ If $b = 3a-\sqrt{6a^2+1}$, then we get $a-\sqrt{6a^2+1} > 0,$ which is a contradiction. Hence we must have $b = 3a + \sqrt{6a^2+1}$ and so $c=b-a = 2a+\sqrt{6a^2+1}$.
At this point, we can describe all solutions, but we still haven't determined exactly when $6a^2+1$ is a perfect square. We've already shown that all solutions to (0) are of the form $(a,b,c) = (a,3a+\sqrt{6a^2+1}, 2a+\sqrt{6a^2+1})$ where $a>1$. Hence all solutions to (a) are of the form $x=3a^2+ak, y=2a^2+ka,$ where $k=\sqrt{6a^2+1}$ is an integer. Conversely, suppose $x=3a^2+ak, y=2a^2+ka,$ where $k=\sqrt{6a^2+1}, $ is an integer and $a > 1.$ Then $ 2x^2+x = 2(3a^2+ak)^2+3a^2+ak = 2(9a^4+6a^3k + a^2k^2) + 3a^2+ak = 18a^4+k(12a^3+a) +3a^2+2a^2(6a^2+1) = 30a^4+5a^2+(12a^3+a)k.$ Also $y+3y^2 = 2a^2+ka + 3(2a^2+ka)^2 = 2a^2+ka + 3(4a^4+4a^3k + k^2a^2) = 12a^4+(6a^2+1)3a^2 + 2a^2+k(12a^3+a) = 30a^4+5a^2+k(12a^3+a)$.
Hence all solutions are given by the above form.
Now, when exactly is $6a^2 + 1 = x^2$ for some perfect square $x^2$?
The above implies $6a^2 = (x-1)(x+1).$ x must be odd, so the gcd of $x-1$ and $x+1$ is 2. Write $x-1 = 2k$ for some integer k. Then we have $3a^2 = k(2k+2),$ so $a$ is even, say $a=2j$ for some $j\ge 1$. Then $6j^2 = k(k+1).$ k and k+1 are coprime, so the only possibilities are the following: 1) $6 | k$, 2) $ 6|(k+1)$, 3) $2 | k, 3|(k+1)$, 4) $3 | k , 2|(k+1).$ Suppose $1)$ holds so that $j^2 = i(6i+1)$ for some integer $i$. Both $i$ and $6i+1$ must be perfect squares. But this seems a bit circular. Case 2) seems similar. So how can I proceed from here?
A more straightforward approach is to use completing the square to identify the relevant Pell-type equation. To wit:
$2x^2+x=3y^2+y$
$2(x^2+\frac12x)=3(y^3+\frac13y)$
$2(x^2+\frac12x+\frac1{16})=3(y^2+\frac13y+\frac1{36})+\frac1{24}$
where the $1/24$ term is included to make up for the difference between the completion terms. So
$2(x+\frac14)^2-3(y+\frac16)^2=\frac1{24}$
$3(4x+1)^2-2(6y+1)^2=1$
$\color{blue}{(12x+3)^2-6(6y+1)^2=3}$
One solution to the generalized Pell equation in blue is given by $12x_0+3=3,6x_0+1=1$ which is the "trivial" solution $x_0=y_0=0$. From the theory of Pell-type equations we are then sure that a series of additional solutions (not always all of them for all Pell indices, but it does turn out to be all solutions when the index is $6$), for positive $x$ and $y$, will be given by
$(12x_n+3)+(6y_n+1)\sqrt6=(3+\sqrt6)(5+2\sqrt6)^n$
where the coefficients $5$ and $2$ come from the minimal positive solution to $k^2-6l^2=1$. So we develop a recursive relation:
$(12x_{n+1}+3)+(6y_{n+1}+1)\sqrt6=((12x_n+3)+(6y_n+1)\sqrt6)(5+2\sqrt6)$
$x_{n+1}=5x_n+6y_n+2$
$y_{n+1}=4x_n+5y_n+\frac53$
Putting in $x_0=y_0=0$ and iterating then gives the positive solution set
$(x,y)\in\{(2,\frac53),(22,18),(220,\frac{539}3),(2180,1780),...\}$
in which, for positive whole number solutions, we accept only every other one from the above set. This alternation is connected with the fact that the primitive solution to $k^2-6l^2=1$ identified above has $k\equiv-1\not\equiv+1\bmod6$.