Метод Рунге-Кутта 2го порядка

#include <iostream>
#include <vector>

using namespace std;


vector<double> y = { 1, 1, 1 };


vector<double> k1(3, 0), k2(3, 0), k3(3, 0), k4(3, 0);

double e = 0.001;

vector<vector<double>> A = {
    { -1.f / e,     1.f / e,           1.f / e },
    { 0,            1.f - 1.f / e,     0 },
    { 0,            0,                 1.f }
};

double h;
double step = 0.1f;

auto getK1() {
    vector<double> result(3, 0);
    for (int i = 0; i < 3; i++) {
        result[i] = step * (A[i][0] * y[0] + A[i][1] * y[1] + A[i][2] * y[2]);
    }
    return result;
}

auto getK2() {
    vector<double> result(3, 0);
    for (int i = 0; i < 3; i++) {
        result[i] = step * (A[i][0] * (y[0] + k1[0]) + A[i][1] * (y[1] + k1[1]) + A[i][2] * (y[2] + k1[2]));
    }
    return result;
}

auto getY() {
    vector<double> result(3, 0);
    for (int i = 0; i < 3; i++) {
        result[i] = y[i] + 1.f / 2.f * (k1[i] + k2[i]);
    }
    return result;
}

int main() {
    for (h = 0 + step; h < 1; h += step) {
        k1 = getK1();
        k2 = getK2();
        y = getY();
        cout << h << ";" << (double)y[0] << ";" << (double)y[1] << ";" << (double)y[2] << endl << endl;
    }
    system("pause");
    return 0;
}