#include <stdio.h>
#include <math.h>
#include <stdlib.h>
 
#define EPSILON 1e-6
#define MAX_ITER 100
 
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;
}
 
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;
    int iter = 0;
 
    printf("Initial guess: x = %.5f, y = %.5f\n", x, y);
 
    while (iter < MAX_ITER) {

        double F1 = f1(x, y);
        double F2 = f2(x, y);

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

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

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

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

  
        printf("iter %2d: x=%.10f  y=%.10f  f1=%.10f  f2=%.10f  dx=%.10f  dy=%.10f  x_next=%.10f  y_next=%.10f\n",
       iter, x, y, F1, F2, dx, dy, x + dx, y + dy);


        x += dx;
        y += dy;
        iter++;
    }
 
    printf("Did not converge within %d iterations.\n", MAX_ITER);
    return 1;
}