I'm trying to solve a system of 20 coupled differential equations through computing. These equations exhibit strong oscillatory solutions over millions of years. I have been reading around, however I can't find a general consensus on which integration methods to use for what, in particular I would like one that has a nice balance between computation time and accuracy.
So far I have been using a dynamic mix of (nonstiff) Adams method/ (stiff) BDF, this is very slow in Python with often inaccurate results over larger timescales.
Would Runge Kutta 4 be a better method for what I am trying to do? How can I generally pick a good method for what I am trying to do?