Метод Ньютона решение СНУ

int n = 2; // Кол-во переменных и функций

double epsilon1, epsilon2;

/**
 * Норма вектора
 */
double norm(std::vector<double> a) {
    double sum = 0;
    for (auto i : a)
        sum += i * i;
    return sqrt(sum);
}

/**
 * Проверка условий, при которых нужно выходить
 * 1. По шагу (т.е. если шаг стал меньше eps1)
 * 2. По норме
 */
int exit_condition(std::vector<double> F, double beta) {
    if (beta < epsilon1)
        return 0;
    if (norm(F) < epsilon2)
        return 0;
    return 1;
}

/**
 * Заданные функции (их кол-во должно быть равно n)
 * X[0] - X[1]        <=>   Функция y = x
 * X[0] + X[1] - 2    <=>   y = -x + 2
 * X[0] - 1           <=>   y = 1
 * Если все функции не пересекаются в одной точке, то одна из функций всегда будет возвращать большое значение в невязке
 * и т.к. эпсилон очень маленький, то решение никогда не сойдется
 * Поэтому нужно брать функции, которые пересекаются в одной точке
 */
double function(std::vector<double> X, int i) {
    if (i == 0) return (X[0] - X[1] * X[1]);
    //if (i == 0) return (X[0] - X[1]);
    if (i == 1) return (X[0] + X[1] - 2);
    if (i == 2) return (X[0] - 1);
}

/**
 * Решаем слау вида  A * DeltaX = F
 * X нужен для вычисления A (матрица Якоби)
 *
 * Составляем матрицу якоби (матрица производных)
 */
void Matrix_Yacoby(std::vector<std::vector<double>> &A, std::vector<double> X, std::vector<double> F, std::vector<double> &DeltaX) {
    double delta = 0.000001;

    // Формула производной f'(x) = ( f(x + Δx) - f(x) ) / Δx
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            X[j] += delta;
            A[i][j] = function(X, i);
            X[j] -= delta;
            A[i][j] -= function(X, i);
            A[i][j] /= delta;
        }
    }
    slove(A, DeltaX, F, n);
}

void input(int &max) {
    std::ifstream A("Text.txt");
    A >> max;
    A >> epsilon1;
    A >> epsilon2;
    A.close();
}

int main() {
    std::vector<std::vector<double>> A(n, std::vector<double>(n, 0));
    std::vector<double> oldF(n), F(n), X(n, 50.0), DeltaX(n, 0.0);

    int amount_iter = 0;
    int max_amount_iter;
    double beta = 1.0;
    input(max_amount_iter);
    for (int i = 0; i < n; i++)
        F[i] = function(X, i);

    while (exit_condition(F, beta) && (amount_iter < max_amount_iter)) {
        beta = 1.0;// начальное значение бета всегда 1
        Matrix_Yacoby(A, X, F, DeltaX); // Считаем вектор DeltaX, т.е. вектор, куда нам нужно идти
        oldF.swap(F);

        for (int i = 0; i < n; i++)
            X[i] -= beta * DeltaX[i];// вычисляем новое приближение Х

        for (int i = 0; i < n; i++)
            F[i] = function(X, i);
        //int k = 0;
        while ((norm(oldF) < norm(F))) {// если норма функции от нового Х больше, чем от старого
            beta /= 2.0;// то уменьшаем бету в 2 раза, вычисляем Х и функцию от него заново
            for (int i = 0; i < n; i++)
                X[i] += beta * DeltaX[i];
            for (int i = 0; i < n; i++)
                F[i] = function(X, i);
        }
        for (auto i : X)
            std::cout << i << " ";
        std::cout << std::endl;

        amount_iter++;
    }

    std::cout << "Amount iteration " << amount_iter << std::endl;
    std::cout << "beta " << beta << std::endl;
    system("pause");
    return 0;
}