Perturbed phase oscillator integration
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.