Mathematica for solving coupled ordinary differential equation

Probably many know that Wolfram Mathematica is a great tool. I knew it as well, but now I'm actually observing first hand how powerful it really is. It is like with chilli pepper - everyone knows that it is hot, but you unless you taste it you won't really understand it. My work involves solving and manipulating many ordinary differential equations (ODE) which quite often are coupled. Writing basic script in Python to do that isn't hard. For simple cases one can use SciPy's build-in function ode from class integrate (documentation). It isn't very fast, and writing everything takes some time, but it was still faster than solution I saw for Matlab or C++. However, the winner, in my opinion, is Mathematica with it's simple structure. I don't have it often, but looking at code for generating and solving $n$ coupled ODEs and seeing how fast it's performed makes me simply happy. Really simple!

Below is code to generate $n$ coupled ODEs and parameters for them. Note: Subscript[x, i] is only syntax and could be exchanged with $x_i$, but for copying purposes I left it in long form.

n = 7;
Do[Do[Subscript[k, i, j] = RandomReal[{0, 1}], {i, 1, n}], {j, 1, n}];
Do[Subscript[w, i] = RandomReal[{2 n, 4 n}], {i, 1, n}];
Do[Subscript[r0, i] = RandomReal[{1, 3}], {i, 1, n}];
Do[Subscript[p, i] = RandomReal[{1, 3}], {i, 1, n}];

eqns = Table[{Subscript[$\phi$, i]'[t] == Subscript[w, i] + Sum[Subscript[k, i, j] Sin[Subscript[$\phi$, j][t] - Subscript[$\phi$, i][t]], {j, 1, n}], Subscript[$\phi$, i][0] == Subscript[p, i]}, {i, 1, n}];

vars = Table[Subscript[$\phi$, i][t], {i, 1, n}];

sol = NDSolve[eqns, vars, {t, 0, 2 Pi}];

Whole magic is in function NDSolve (documentation), which solves numerically equations eqns for variables vars in provided range. Now all what's left is to plot results. Again, very simple.

This code generates plots and display them in column. Output graphs below.

FreqTable = Table[Plot[Evaluate[D[Subscript[$\phi$, i][t] /. sol, t]], {t, 0, 2 Pi}, GridLines -> {{}, {Subscript[w, i]}}, AxesLabel -> {Time[s], InstFreq [1/s]}], {i, 1, n}];
GraphicsColumn[FreqTable, ImageSize -> Full]