import numpy as np

def f(x):
    # Nonlinear system
    x_val, y_val = x[0], x[1]
    f1 = 4 * x_val - 2 * x_val * y_val
    f2 = 2 * y_val + x_val * y_val - 2 * y_val**2
    return np.array([f1, f2])

def J(x):
    # Jacobian of the system
    x_val, y_val = x[0], x[1]
    j11 = 4 - 2 * y_val
    j12 = -2 * x_val
    j21 = y_val
    j22 = 2 + x_val - 4 * y_val
    return np.array([[j11, j12],
                     [j21, j22]])

def newton_method(initial, iterations):
    x = np.array(initial, dtype=float)
    history = [x.copy()]
    for _ in range(iterations):
        F = f(x)
        Jx = J(x)
        # Add small regularization to avoid singular matrix
        reg = 1e-10
        Jx += reg * np.eye(Jx.shape[0])
        # Solve J * dx = F for dx, then update
        dx = np.linalg.solve(Jx, F)
        x = x - dx
        history.append(x.copy())
    return history

# Run Newton's method for 20 iterations and print the last 10 iterations
if __name__ == '__main__':
    iterations = 20
    history = newton_method([1, 1], iterations)
    for i, sol in enumerate(history[-10:], start=len(history)-10):
        print(f"Iteration {i}: x = {sol[0]}, y = {sol[1]}")