Solve the system in $mathbb{R}$

Using your approach, consider the function $$f=k^5+2k^4+193k^3-1044k^2+2160k-1728$$ which has only one real zero close to $k=2$.

Rewrite $f$ as a Taylor series (!!) around $k=2$ to get $$f=24+444 (k-2)+242 (k-2)^2+249 (k-2)^3+12 (k-2)^4+(k-2)^5+Oleft((k-2)^6right)$$ and use series reversion to get $$k=2+t-frac{121 }{222}t^2+frac{1643 }{49284}t^3+frac{3784043 }{5470524}t^4-frac{1497968351}{1214456328} t^5+Oleft(t^6right)$$ where $t=frac{f-24}{444}$.

Make $f=0$ (that is to say $t=-frac 2 {37}$) to get $$k=frac{20468028448605260}{10526900923401237}approx 1.944354620$$ while the exact solution is $1.944354635$.

begin{align} x+y+frac{x^2}{y^2}=7 tag{1}label{1} ,\ (x-y)frac{x^2}{y^2}=12 tag{2}label{2} . end{align}

Let $tfrac xy=t$. Then we can transform eqref{1}-eqref{2} into linear system of two variables $x,y$ in terms of parameter $t$:

begin{align} x+y &=7-t^2 tag{3}label{3} ,\ x-y &= frac{12}{t^2} tag{4}label{4} ,\ x&= -frac{t^4-7t^2-12}{2t^2} tag{5}label{5} ,\ y&= -frac{t^4-7t^2+12}{2t^2} tag{6}label{6} . end{align}

Dividing eqref{5} by eqref{6}, we get

begin{align} t &= frac{t^4-7t^2-12}{t^4-7t^2+12} tag{7}label{7} , end{align}

which is equivalent to irreducible quintic equation begin{align} t^5-t^4-7t^3+7t^2+12t+12 &=0 tag{8}label{8} end{align}

with one real root begin{align} tauapprox -2.4842933862 tag{9}label{9} , end{align}

hence the original system has one real solution begin{align} x&= -frac{tau^4-tau^2-12}{2tau^2} tag{10}label{10} approx 1.386320 ,\ y&= -frac{tau^4-7tau^2+12}{2tau^2} tag{11}label{11} approx -0.558034 . end{align}

Method I: Nonlinear Newton's Method for a system of equations.

Using Solving a set of equations with Newton-Raphson, the regular Newton-Raphson method is initialized with a starting point $x_0$ and then you iterate $tag 1x_{n+1}=x_n-dfrac{f(x_n)}{f'(x_n)}$

In higher dimensions, there is an exact analog. We define:

$$Fleft(begin{bmatrix}x_1\x_2end{bmatrix}right) = begin{bmatrix}f_1(x_1,x_2) \ f_2(x_1,x_2) end{bmatrix} =begin{bmatrix}dfrac{x_1^2}{x_2^2}+x_1+x_2-7,dfrac{x_1^2 (x_1-x_2)}{x_2^2}-12end{bmatrix}$$

The derivative of this system is the $2x2$ Jacobian given by:

$$J(x_1,x_2) = begin{bmatrix} dfrac{2 x_1}{x_2^2}+1 & 1-dfrac{2 x_1^2}{x_2^3} \ dfrac{x_1^2}{x_2^2}+dfrac{2 x_1 (x_1-x_2)}{x_2^2} & -dfrac{x_1^2}{x_2^2}-dfrac{2 x_1^2 (x_1-x_2)}{x_2^3} end{bmatrix}$$

The function $G$ is defined as:

$$G(x) = x - J(x)^{-1}F(x)$$

and the functional Newton-Raphson method for nonlinear systems is given by the iteration procedure that evolves from selecting an initial $x^{(0)}$ and generating for $k ge 1$ (compare this to $(1)$),

$$x^{(k)} = G(x^{(k-1)}) = x^{(k-1)} - J(x^{(k-1)})^{-1}F(x^{(k-1)}).$$

We can write this as:

$$begin{bmatrix}x_1^{(k)}\x_2^{(k)}end{bmatrix} = begin{bmatrix}x_1^{(k-1)}\x_2^{(k-1)}end{bmatrix} + begin{bmatrix}y_1^{(k-1)}\y_2^{(k-1)}end{bmatrix}$$

where:

$$begin{bmatrix}y_1^{(k-1)}\y_2^{(k-1)}end{bmatrix}= -left(Jleft(x_1^{(k-1)},x_2^{(k-1)}right)right)^{-1}Fleft(x_1^{(k-1)},x_2^{(k-1)}right)$$

Here is a starting vector:

$$x^{(0)} = begin{bmatrix}x_1^{(0)}\x_2^{(0)}end{bmatrix} = begin{bmatrix}1.0\-1.0end{bmatrix}$$

You should end up with a solution that looks something like (you can do the iteration formula steps as it converges very quickly):

$$begin{bmatrix}x_1\x_2end{bmatrix} = begin{bmatrix} 1.38632 \ -0.558034 end{bmatrix}$$

Method II: Use a Groebner Basis to eliminate a variable in the equations.

  GroebnerBasis[{x + y + x^2/y^2 - 7, (x - y) x^2/y^2 - 12}, {x,  y}, {y}]

Result of that is $$x^5-9 x^4+96 x^3-396 x^2+1008 x-864 = 0$$

We could have also eliminated $x$ and solved

$$y^5-11 y^4+79 y^3-33 y^2+84 y+72 = 0$$

Using numerical methods on either of those, gives the same result as above.

Method III: The approach you outlined. The system is given as

$$begin{align} x+y+dfrac{x^2}{y^2}&=7 \ dfrac{(x-y)x^2}{y^2}&=12 end{align}$$

Using the second equation, we solve as $frac{x^2}{y^2} = frac{12}{x-y}$ and substitute into the first

$$(x-y)(x+y) + 12 = 7(x-y)$$

We now choose $$x = y+k$$ and substitute and arrive at

$$y = dfrac{-k^2 + 7k-12}{2k} implies x = dfrac{-k^2 + 7k-12}{2k} + k$$

Substituting these into the first equation, we have

$$dfrac{k^5+2 k^4+193 k^3-1044 k^2+2160 k-1728}{k left(k^2-7 k+12right)^2} = 0$$

Using Newton's Method on the numerator yields

$$k = 1.9443546350028582064014095438170755609716945099813650490087$$

Substituting this $k$ into the values for $x$ and $y$ yields the same results as shown above.

Notice that in all three cases, we had to resort to numerical methods, the first using Nonlinear Newton's and the next two using Newton's Method because of an irreducible quintic.