#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.001f;

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.0001f;

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] + 0.5f * (k1[i] + k2[i]);
	}
	return result;
}

int main() {
	double d = 1 / step;
	d = d / 10 + 0.1;
	int n = 1, p = d;
	for (h = 0 + step; h < 1; h += step, n++) {
		k1 = getK1();
		k2 = getK2();
		y = getY();
		if ((n % p) == 0)
			cout << h << "; " <<(double)y[0] << "; " <<(double)y[1] << "; " <<(double)y[2] << endl << endl;
	}
	system("pause");
	return 0;
}