#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define EPSILON 1e-6   // Convergence tolerance
#define MAX_ITER 100   // Maximum iterations

// Define the nonlinear equations
double f1(double x, double y) {
    return x * x + y * y - 4; 
}

double f2(double x, double y) {
    return x * x * x * x + x - y*y + 1; // Hyperbola: xy = 1
}

// Partial derivatives for Jacobian matrix
double df1_dx(double x, double y) { return 2 * x; }
double df1_dy(double x, double y) { return 2 * y; }
double df2_dx(double x, double y) { return 4 * x * x*x + 1; }
double df2_dy(double x, double y) { return -2* y; }

int main() {
    double x = 1, y = 3; // Initial guess
    int iter = 0;

    printf("Initial guess: x = %.6f, y = %.6f\n", x, y);

    while (iter < MAX_ITER) {
        // Function values
        double F1 = f1(x, y);
        double F2 = f2(x, y);

        // Check convergence
        if (fabs(F1) < EPSILON && fabs(F2) < EPSILON) {
            printf("Converged after %d iterations.\n", iter);
            printf("Solution: x = %.6f, y = %.6f\n", x, y);
            return 0;
        }

        // Jacobian matrix
        double J11 = df1_dx(x, y);
        double J12 = df1_dy(x, y);
        double J21 = df2_dx(x, y);
        double J22 = df2_dy(x, y);

        // Determinant of Jacobian
        double det = J11 * J22 - J12 * J21;
        if (fabs(det) < 1e-12) {
            printf("Jacobian determinant too small. Method fails.\n");
            return 1;
        }

        // Inverse Jacobian * F
        double dx = (-F1 * J22 + F2 * J12) / det;
        double dy = (-J11 * F2 + J21 * F1) / det;

        // Update variables
        x += dx;
        y += dy;

        iter++;
    }

    printf("Did not converge within %d iterations.\n", MAX_ITER);
    return 1;
}