Как улучшить скорость odeint в Python?

Я использую Python и odeint из пакета scipy для решения большого числа ( ~ 10e6) связанных ОДУ. Систему уравнений можно сформулировать как сумму некоторых умножений матриц, и я использую для этого numpy с поддержкой blas. Моя проблема в том, что это занимает очень много времени. Когда я профилирую код, я вижу, что большинство времени уходит на то, чтобы делать что-то еще, кроме оценки правых. Это пять наиболее трудоемких вызовов от профилировщика:

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
5       1547.915  309.583 1588.170  317.634 {scipy.integrate._odepack.odeint}
60597   11.535    0.000   23.751    0.000   terms3D.py:5(two_body_evolution)
121194  11.242    0.000   11.242    0.000   {numpy.core._dotblas.dot}
60597   10.145    0.000   15.460    0.000   generator.py:13(Gs2)
121203   3.615    0.000   3.615     0.000   {method 'repeat' of 'numpy.ndarray' objects}

Правая часть состоит в основном из two_body_evolution и Gs2. Этот профиль e для ~ 7000 связанных ODE, и вот то же самое для ~ 4000:

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
5        259.427  51.885   273.316  54.663 {scipy.integrate._odepack.odeint}
30832    3.809    0.000    7.864    0.000  terms3D.py:5(two_body_evolution)
61664    3.650    0.000    3.650    0.000  {numpy.core._dotblas.dot}
30832    3.464    0.000    5.637    0.000  generator.py:13(Gs2)
61673    1.280    0.000    1.280    0.000  {method 'repeat' of 'numpy.ndarray' objects}

Итак, моя главная проблема здесь в том, что «скрытое» время в odeint ужасно масштабируется с количеством уравнений. У вас есть идеи, почему это так и как улучшить производительность?

Спасибо за ваше время

Оскар Окерлунд