如果你能用Matlab的ode15s解决你的问题，你应该也能用scipy的vode求解器来解决它。为了模拟ode15s，我使用了以下设置：

ode15s = scipy.integrate.ode(f)
ode15s.set_integrator('vode', method='bdf', order=15, nsteps=3000)
ode15s.set_initial_value(u0, t0)

然后你就可以很高兴地用ode15s.integrate(t_final)解决你的问题了。它在一个僵硬的问题上应该能很好地工作。
(See也称为Link）

Python可以调用C。ODEPACK中的行业标准是LSODE。它是公共领域的。你可以下载C version。这些解算器非常棘手，所以最好使用一些经过良好测试的代码。
补充：确保你真的有一个刚性系统，也就是说，如果速率（特征值）相差超过2或3个数量级。此外，如果系统是刚性的，但你只是在寻找一个稳态解，这些解算器给予你的选择，解决一些方程代数。否则，一个好的龙格库塔解算器，如DVERK将是一个很好的，更简单的解决方案。
在此添加，因为它不适合放在注解中：这来自DLSODE标题文档：

C     T     :INOUT  Value of the independent variable.  On return it
C                   will be the current value of t (normally TOUT).
C
C     TOUT  :IN     Next point where output is desired (.NE. T).

当然，米氏动力学是非线性的，但艾特肯加速度也适用。（如果您需要简短的说明，首先考虑Y是标量的简单情况。运行系统得到3个Y（T）点。通过它们拟合指数曲线（简单代数）。然后将Y设为渐近线并重复。现在将Y推广为向量。假设3个点在一个平面上-如果它们不在平面上也可以。）此外，除非你有一个强制功能（像一个恒定的静脉滴注），MM消除将衰减消失，系统将接近线性。希望这有帮助。

PyDSTool封装了Radau解算器，它是一个优秀的隐式刚性积分器。它比odeint有更多的设置，但比PyGSL少很多。最大的好处是你的RHS函数被指定为字符串（通常，尽管你可以使用符号操作来构建系统），并被转换为C，所以没有缓慢的python回调，整个过程将非常快。

我目前正在研究一点ODE和它的求解器，所以你的问题对我来说很有趣...
从我所听到的和读到的，对于僵硬的问题，正确的方法是选择一个隐式方法作为阶跃函数（如果我错了请纠正我，我还在学习ODE求解器的奥秘）。我不能引用你我在哪里读到的，因为我不记得了，但这里有一个来自gsl-help的thread，在那里问了一个类似的问题。
因此，简而言之，bsimp方法似乎值得一试，尽管它需要雅可比函数。如果你不能计算雅可比函数，我将尝试rk2imp、rk4imp或任何齿轮方法。