I try to simplify a polynomial to the form: $y^2=4x^3-g_{2}x-g_{3}$, which is the elliptic curves. And the polynomial is $xy(x+y)+(x+y)+(x+y)^2=13xy$.
I try to let the $u=x+y$ and $v=x-y$, then I get the $u^3-uv^2+4u-9u^2+13v^2=0$. But how to get ahead?
Most working mathematicians would use a computer program for this computation. Magma works well, and is available at no charge in all US universities. Below is the computation in Magma. Note that Magma requires the use of projective coordinates. To convert from projective to affine, just set $Z=1$.
> P2<X,Y,Z> := ProjectiveSpace(Rationals(),2);
> C := Curve(P2, X*Y*(X+Y) + (X+Y)*Z^2 + (X+Y)^2*Z - 13*X*Y*Z);
> pt := C![0,0,1];
> EllipticCurve(C,pt);
Elliptic Curve defined by y^2 = x^3 + 121*x^2 + 1144*x + 2704 over Rational Field
Mapping from: CrvPln: C to Elliptic Curve defined by y^2 = x^3 + 121*x^2 + 1144*x + 2704 over Rational Field
with equations :
52*X*Z
-52*X^2 - 52*X*Y + 572*X*Z - 104*Y*Z - 104*Z^2
-X^2 - X*Y
What all of the above means is that the elliptic curve $y^2 = x^3 + 121x^2 + 1144x + 2704$ under the substitutions \begin{align*} x &= \frac{52 X}{-X^2 - X Y} = \frac{52}{-X-Y}\\ y &= \frac{-52X^2 - 52XY + 572X - 104Y - 104}{-X^2 - X Y} \end{align*} yields, after some simplification, $$ \frac{10816 (X + Y + 1)(XY(X+Y) + (X+Y) + (X+Y)^2 - 13XY)}{X^2(X+Y)^3} = 0, $$ which is equivalent to your original curve (look in the numerator of the fraction above).
The elliptic curve $y^2 = x^3 + 121x^2 + 1144x + 2704$ isn't quite in the form that you want. Converting it into the form that you want is an easy exercise.
Sage (an open-source alternative to Magma) also works, but less well. It returns enormous equations:
sage: R.<x,y> = QQ[]
sage: cubic = x*y*(x+y) + (x+y) + (x+y)^2 - 13*x*y
sage: EllipticCurve(WeierstrassForm(cubic))
Elliptic Curve defined by y^2 = x^3 + (-11209/48)*x + 1185157/864 over Multivariate Polynomial Ring in x, y over Rational Field
sage: WeierstrassForm(cubic,transformation=True)
(18421/8748*x^6*y^6 - 2873/1458*x^6*y^5 - 2873/1458*x^5*y^6 + 1692157/69984*x^6*y^4 + 1124317/34992*x^5*y^5 + 1692157/69984*x^4*y^6 + 14351/1296*x^6*y^3 - 5844665/17496*x^5*y^4 - 5844665/17496*x^4*y^5 + 14351/1296*x^3*y^6 + 1692157/69984*x^6*y^2 - 5844665/17496*x^5*y^3 + 150149899/69984*x^4*y^4 - 5844665/17496*x^3*y^5 + 1692157/69984*x^2*y^6 - 2873/1458*x^6*y + 1124317/34992*x^5*y^2 - 5844665/17496*x^4*y^3 - 5844665/17496*x^3*y^4 + 1124317/34992*x^2*y^5 - 2873/1458*x*y^6 + 18421/8748*x^6 - 2873/1458*x^5*y + 1692157/69984*x^4*y^2 + 14351/1296*x^3*y^3 + 1692157/69984*x^2*y^4 - 2873/1458*x*y^5 + 18421/8748*y^6,
15379/6561*x^9*y^8 - 15379/6561*x^8*y^9 - 15379/729*x^9*y^7 + 15379/729*x^7*y^9 + 1491763/19683*x^9*y^6 - 61516/6561*x^8*y^7 + 61516/6561*x^7*y^8 - 1491763/19683*x^6*y^9 - 353717/6561*x^9*y^5 - 61516/2187*x^8*y^6 + 61516/2187*x^6*y^8 + 353717/6561*x^5*y^9 + 353717/6561*x^9*y^4 + 76895/6561*x^7*y^6 - 76895/6561*x^6*y^7 - 353717/6561*x^4*y^9 - 1491763/19683*x^9*y^3 + 61516/2187*x^8*y^4 - 76895/6561*x^7*y^5 + 76895/6561*x^5*y^7 - 61516/2187*x^4*y^8 + 1491763/19683*x^3*y^9 + 15379/729*x^9*y^2 + 61516/6561*x^8*y^3 + 76895/6561*x^6*y^5 - 76895/6561*x^5*y^6 - 61516/6561*x^3*y^8 - 15379/729*x^2*y^9 - 15379/6561*x^9*y - 61516/6561*x^7*y^3 - 61516/2187*x^6*y^4 + 61516/2187*x^4*y^6 + 61516/6561*x^3*y^7 + 15379/6561*x*y^9 + 15379/6561*x^8*y - 15379/729*x^7*y^2 + 1491763/19683*x^6*y^3 - 353717/6561*x^5*y^4 + 353717/6561*x^4*y^5 - 1491763/19683*x^3*y^6 + 15379/729*x^2*y^7 - 15379/6561*x*y^8,
-13/27*x^3*y^3 + 121/108*x^3*y^2 + 121/108*x^2*y^3 + 121/108*x^3*y - 1175/108*x^2*y^2 + 121/108*x*y^3 - 13/27*x^3 + 121/108*x^2*y + 121/108*x*y^2 - 13/27*y^3)
Are you certain that you have stated the problem correctly? As currently stated
$$ xy(x+y)+(x+y)+(x+y)^2=13xy $$
can be rewritten
\begin{eqnarray} x^2y+xy^2+x+y+x^2+2xy+y^2&=&13xy\\ y^2(x+1)+y(x^2-11x+1)+x&=&0\\ y(x^2-11x+1)&=&-[y^2(x+1)+x]\\ y^2(x^2-11x+1)^2&=&[y^2(x+1)+x]^2 \end{eqnarray} So if there exist constants $g_3,\,g_2$ such that $$y^2=4x^3-g_{2}x-g_{3}$$
then
$$ (4x^3-g_{2}x-g_{3})(x^2-11x+1)^2=[(4x^3-g_{2}x-g_{3})(x+1)+x]^2 $$
This will resolve in an eighth degree polynomial equation of the form
$$ 16x^8-4x^7+P_6(x,g_2,g_3)=0 $$
Where $P_6$ is a degree six polynomial whose cofficients are function of $g_3,\,g_2$. But the coefficients of $x^8$ and $x^7$ are independent of the values of $g_3,\,g_2$.
So no constant values of $g_3,\,g_2$ can result in a solution for $x\ne0$.
Thus there cannot be values of $g_3,\,g_2$ for which $y^2=4x^3-g_{2}x-g_{3}$.
You made a good start but there was a complication you did not anticipate. The homogeneous version of your equation is: $\, 0 = -W X Y Z + (X+Y)(X+Z)(Y+Z), \,$ where $\, W=13 \,$ is a constant. Now substitute $\, X = 1 + c_3 x + \sqrt{c_1} y, $ $\, Y = 1 + c_3 x - \sqrt{c_1} y, $ $\, Z = 2 x + c_2, \,$ where $\,c_1,c_2,c_3 \,$ depend on $\,W.\,$ After the substitutions, we eliminate the $\, x y^2 \,$ term with $\, c_3 = W. \,$ We eliminate the $\,x^2\,$ term with $\, c_2 = (2W^2 +16W + 8)/(W^3 - 4W^2 - 8W). \,$ Now the coefficients of $\,y^2\,$ and of $\,x^3\,$ needs adjustment to get the final form. Let $\, c_1 = (-W^3 + 4W^2 + 8W)/12 .\,$ The equation now is $\, 0 = 2W(1+W)( -y^2 + 4x^3 - g_2x - g_3), \,$ where $$\, g_2 = \frac{12(W^4 - 8W^3 + 16W + 16)}{W^2(W^2 - 4W - 8)^2}, \quad g_3 = \frac{8(W^4 - 8W^3 - 8W - 8)}{W^3(W^2 - 4W - 8)^2}. $$ In our case of $\, W=13 \,$ these invariants become $\, g_2 = 134508/2007889, \, g_3 = 86984/26102557. \,$ There is a $2$-torsion point $\, (-1/13,0). \,$ There is a generator point $\, (-97/1417, 168 \sqrt{-3}/1417^{3/2}). \,$ After noticing the $1417$ in the denominators, we can scale $\,x,y\,$ to simplify the equation. After scaling it is $\, 0 = -y^2 + 4x^3 - 134508x - 9481256, \,$ with generator point $\, (-97, 168 \sqrt{-3}). \,$
This is the elliptic curve with LMFDB label 8190.bw4 which has no rational generator.
