Runge-Kutta fourth order integration I have written a simple code for Runge-Kutta fourth order integration to solve a system of ordinary differential equations and parallelized it using OpenMP. I don't know if it is the best we can do for maximum performance of the code with little effort.

I need all values of to be returned, so I kept values in all steps.
I also create threads in each time step and paralleled in position, i.e. 
pragma opm parallel is inside the loop over time.

Here is my try: (link to gitlab repository)

//RHS for a system of equations
void xprsys(const int n,const vector&lt;double&gt;&amp; x, vector&lt;double&gt;&amp; f)
{
    /*
     * n : number of equations
     * x : value in each time step
     * f : RHS of equtions dx/dt = f 
     */
    double sum1=0;
    #pragma omp parallel for reduction(+:sum1)
    for (int i=0; i&lt;n; i++){
        sum1 = 0;
        for(int j=0; j&lt;n; j++)
            sum1 += sin(x[j]-x[i]);
        f[i] = M_PI + 2.0 * sum1;
    }
}

void SolveRK4(const int n, double h,vector&lt;double&gt; x,
     vector&lt;vector&lt;double&gt;&gt;&amp; x_vec, 
          int nstep, vector&lt;double&gt;&amp; times)
{
    /*
     * times : vector contains the time 0 : t_final step dt
     * x_vec : [nstep by N] 2 Dimensional vector
     * dim1  : defined as typedef vector&lt;double&gt; dim1
     */

    times[0] = 0.0;
    dim1 y(n);
    dim1 f1(n),f2(n),f3(n),f4(n);
    double half_h = 0.5 * h;
    double h_sixth = h/6.0;
    // x_vec[nstepxN]

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

    for (int k=1; k&lt;nstep; k++){
        times[k] = k*h;
        xprsys(n,x,f1);

        #pragma omp parallel for 
        for(int i=0; i&lt;n; i++)
            y[i] = x[i] + half_h * f1[i];

        #pragma omp master
        xprsys(n,y,f2);
        #pragma omp barrier

        #pragma omp for
        for(int i=0; i&lt;n; i++)
            y[i] = x[i] + half_h * f2[i];

        #pragma omp master
        xprsys(n,y,f3);
        #pragma omp barrier

        #pragma omp for
        for(int i=0; i&lt;n; i++)
            y[i] = x[i] + h * f3[i];

        #pragma omp master
        xprsys(n,y,f4);
        #pragma omp barrier

        #pragma omp for 
        for(int j=0; j&lt;n; j++) {
            x[j] = x[j] + h_sixth * (f1[j] + f4[j] + 2.0 * (f2[j] + f3[j]));
            x_vec[k][j] = x[j];
        }
    }
}
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.