Any way to solve a system of coupled differential equations in python?

For the numerical solution of ODEs with scipy, see scipy.integrate.solve_ivp, scipy.integrate.odeint or scipy.integrate.ode.

Some examples are given in the SciPy Cookbook (scroll down to the section on "Ordinary Differential Equations").

In addition to SciPy methods odeint and ode that were already mentioned, it now has solve_ivp which is newer and often more convenient. A complete example, encoding [v11, v22, v12] as an array v:

from scipy.integrate import solve_ivpdef rhs(s, v):     return [-12*v[2]**2, 12*v[2]**2, 6*v[0]*v[2] - 6*v[2]*v[1] - 36*v[2]]res = solve_ivp(rhs, (0, 0.1), [2, 3, 4])

This solves the system on the interval (0, 0.1) with initial value [2, 3, 4]. The result has independent variable (s in your notation) as res.t:

array([ 0.        ,  0.01410735,  0.03114023,  0.04650042,  0.06204205,        0.07758368,  0.0931253 ,  0.1       ])

These values were chosen automatically. One can provide t_eval to have the solution evaluated at desired points: for example, t_eval=np.linspace(0, 0.1).

The dependent variable (the function we are looking for) is in res.y:

array([[ 2.        ,  0.54560138,  0.2400736 ,  0.20555144,  0.2006393 ,         0.19995753,  0.1998629 ,  0.1998538 ],       [ 3.        ,  4.45439862,  4.7599264 ,  4.79444856,  4.7993607 ,         4.80004247,  4.8001371 ,  4.8001462 ],       [ 4.        ,  1.89500744,  0.65818761,  0.24868116,  0.09268216,         0.0345318 ,  0.01286543,  0.00830872]])

With Matplotlib, this solution is plotted as plt.plot(res.t, res.y.T) (the plot would be smoother if I provided t_eval as mentioned).

Finally, if the system involved equations of order higher than 1, one would need to use reduction to a 1st order system.