As with integration, SciPy has some extremely accurate general-purpose solvers for systems of ordinary differential equations of first order:

For real-valued functions, we have basically two flavors: ode (with options passed with the set_integrator method) and odeint (simpler interface). The syntax of ode is as follows:

ode(f,jac=None)

The first parameter, f, is the function to be integrated, and the second parameter, jac, refers to the matrix of partial derivatives with respect to the dependent variables (the Jacobian). This creates an ode object, with different methods to indicate the algorithm to solve the system (set_integrator), the initial conditions (set_initial_value), and different parameters to be sent to the function or its Jacobian.

The options for integration algorithm are 'vode' for real-valued variable coefficient ODE solver, with fixed-leading-coefficient implementation (it provides Adam's method for non-stiff problems and BDF for stiff); 'zvode...