Решение системы обыкновенных дифференциальных уравнений (не менее 4-х)

#include <iostream> 
#include <functional>
 #include <vector> 
#include <iomanip> 
using namespace std; 
double dx1dt(const double t, const double x1, const double x2, const double x3) {  return x1 - x2 + x3; } 
 double dx2dt(const double t, const double x1, const double x2, const double x3) { return x1 + x2 - x3; } 
 double dx3dt(const double t, const double x1, const double x2, const double x3) { return 2*x1 - 2*x2; } 
 void rungeKutta(double x, double y, double z, double t0, double max_t, double h) { 
double k1[3], k2[3], k3[3], k4[3]; 
double t = t0; 
vector<function<double(double, double, double, double)>> vfunc;
vfunc.push_back(dx1dt); 
vfunc.push_back(dx2dt); 
vfunc.push_back(dx3dt); 
 size_t max_iter = (max_t - t0) / h; 
for (int i = 0; i < max_iter; i++) { 
 for (int j(0); j < 3; ++j) { k1[j] = h*vfunc[j](t, x, y, z);  } 
 for (int j(0); j < 3; ++j) { k2[j] = h*vfunc[j](t + 0.5*h, x + 0.5*k1[0], y + 0.5*k1[1], z + 0.5*k1[2]);  } 
 for (int j(0); j < 3; ++j) { k3[j] = h*vfunc[j](t + 0.5*h, x + 0.5*k2[0], y + 0.5*k2[1], z + 0.5*k2[2]);  } 
for (int j(0); j < 3; ++j) { k4[j] = h*vfunc[j](t + h, x + k3[0], y + k3[1], z + k3[2]);  } 
 x += (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.; 
 y += (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.; 
 z += (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) / 6.;
t = t + h; 
} 
} 
 int main() { double x0 = 1; double y0 = 1; double z0 = 1; double beg = 0; double max_time = 1; double step = 0.1;
 rungeKutta(x0, y0, z0, beg, max_time, step);
 }