Perturbed phase oscillator integration I am integrating a system of perturbed phase oscillators. I define the system of equation and also the Jacobian matrix. I have to reshape the one dimensional vector of system state to a 2 dimensional vector then do the matrix production.

Is there any single-threaded way to make the program faster?

using dim1 = std::vector&lt;double&gt;;
using dim2 = std::vector&lt;std::vector&lt;double&gt;&gt;;

dim2 multiply_matrices(dim2 &amp;A, dim2 &amp;B)
{
    int r1 = A.size();
    int c1 = A[0].size();
    int r2 = B.size();
    int c2 = B[0].size();

    assert(c1==r2 &amp;&amp; "bad dimension for matrix multiplication");

    dim2 mult(r1, dim1(c2));

    for (int i = 0; i &lt; r1; ++i) {
        for (int j = 0; j &lt; c2; ++j) {
            mult[i][j] = 0;
        }
    }

    // Multiplying matrix A and B and storing in array mult.
    for (int i = 0; i &lt; r1; ++i) {
        for (int j = 0; j &lt; c2; ++j) {
            for (int k = 0; k &lt; c1; ++k) {
                mult[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return mult;
}

class kuramoto_exe 
{
    vector&lt;vector&lt;double&gt;&gt; &amp;adj;
    dim1 &amp;omega;
    int kind;
    int N;
    double coupling;

public:
    kuramoto_exe (
        int N_, int kind_, double coupling_,
        dim1&amp; omega_, dim2&amp; adj_) : adj(adj_), omega(omega_), 
            kind(kind_), coupling(coupling_), N(N_) { }

        void operator() (const dim1 &amp;X, dim1&amp; dXdt, const double /*t*/)
    {


        dim1 x(X.begin(), X.begin() + N); // initial condition

        // reshaping:every colomn is a vector (reshape order "F")
        dim2 Y(N, dim1(N));
        for (int j=0; j&lt;N; j++)
            for (int i=0; i&lt;N; i++)
                Y[i][j] = X[N*j+i+N];  


        dim2 Jac(N, dim1(N));

        // definition of unperturbed model
        for (int i=0; i&lt;N; i++) {
            double sumj =0.0;
            for (int j=0; j&lt;N; j++) {
                if ((i != j) &amp;&amp; (adj[i][j] &gt; 1e-8))
                        sumj += adj[i][j] * sin(x[j]-x[i]);

                dXdt[i] =  omega[i]+ coupling * sumj;
            }
        }
        // calculation of the Jacobian 
        for (int i=0; i&lt;N; i++) {
            for(int j=0; j&lt;N; j++) {
                if (i!=j) 
                    Jac[i][j] = coupling * adj[i][j] * cos(x[j] - x[i]);
                else {
                    double sumj = 0.0;
                    for (int jj=0; jj&lt;N; jj++) {
                        if (jj!=i)
                            sumj += adj[i][jj] * cos(x[jj] - x[i]);
                    }
                    Jac[i][j] = -1.0 * coupling * sumj;
                }
            }
        }

        dim2 M = multiply_matrices(Jac, Y);

        for (int j=0; j&lt;N; j++)
            for (int i=0; i&lt;N; i++)
                dXdt[N*j+i+N] = M[i][j]; // revese of reshape order "F"
    }
};
First off, be sure you're compiling with optimizations enabled.

In multiply_matrices you don't need to explicitly fill your matrix with zeros, as the vector constructor will do that.

In the kuramoto_exe constructor, the member initializer list has the member variables in the wrong order.  They will be constructed in the order they are declared, so N will be constructed before coupling despite its appearing later in the list.  While not an issue here, it can be if there are dependencies among the member variables, and some compilers will issue a warning when you do this.

In your matrix multiply, accumulate the result to a local variable, then store the result in mult[i][j].  This avoids the potential of multiple index calculations.

When calculating the Jacobian, the diagonal members are the negative of the sum of the horizontal members.  You're doing twice as much work as you need to.  Accumulate the sum as you process the row (skipping the diagonal element) then store the diagonal element when you're done with the row (Jac[i][i] = -rowsum;).

When you're doing the initial reshaping (filling the Y matrix), it is probably better to put the i loop on the outside, so that your writes are sequential.  It is easier for the CPU to handle nonsequential reads then nonsequential writes.

Rather than using a vector of vectors, consider creating a matrix class that stores all elements contiguously (in one vector) and computes the appropriate index.  This avoids the double lookup you currently have.