Runge-Kutta fourth order integration
I would really suggest that you better structure your code. Runge Kutta methods are quite straightforward, so you should have a function that iterates a single step in parallel for every variable. 

Furthermore, you can slightly deviate from the standard Runge Kutta scheme if you directly store \$x + h_if_i\$ in the \$i\$th component of that vector. You only have to adopt the summation at the end to account for that.

I did so here https://github.com/miscco/Neural_Network_TC/blob/master/Iterate_ODE.h

static void add_RK(std::vector&lt;double&gt;&amp; var) {
    var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6;
}


A minimal improvement that I would do is changing this

for (int i=0; i&lt;n; i++)
    x_vec[0][i] = x[i];


To 

x_vec[0] = x;


Where you utilize the underlying copy/move constructor, which might be more efficient depending on the compiler.