先看你给出的测试代码：这个测试里的ODE是dx/dt = 0，解是恒定值0.5*M0，积分区间从t0=0.5*M0到t1=M0，步长dt=0.01*M0，总共有(M0 - 0.5*M0)/0.01*M0 = 50步，所以数组bli有50个元素是完全正常的——这说明在这个简单场景下，网格量级的设置并没有问题，那你原ODE求解失败的原因可能出在其他和数值尺度、求解器适配相关的点上。

下面针对你怀疑的「网格量级」相关方向，结合SciPy的ode求解器给出具体排查建议：

1. 检查原ODE是否为刚性系统 #

你当前使用的dopri5是自适应步长的非刚性求解器，如果你的原ODE是刚性系统（方程中同时存在极快和极慢的变化尺度），非刚性求解器会因为需要极小步长而出现求解超时、数值溢出等问题。这种情况下，建议更换为刚性求解器：

# 方法1：使用lsoda（自动适配刚性/非刚性系统）
r = ode(your_ode_func).set_integrator('lsoda')

# 方法2：使用vode并指定BDF方法（专门处理刚性系统）
r = ode(your_ode_func).set_integrator('vode', method='bdf')

2. 避免强制固定步长，发挥求解器的自适应能力 #

你的测试代码中用r.integrate(r.t+dt)强制固定步长，但dopri5本身是自适应步长求解器，强制固定步长会限制它根据方程变化调整步长的能力，可能在原ODE的复杂区域出错。建议尝试直接积分到终点，让求解器自动选择合适步长：

r.set_initial_value(y0, t0)
if r.successful():
    r.integrate(t1)
    print("最终解：", r.y)

如果需要分步获取结果，可以通过r.t追踪当前积分时间，或调整求解器的nsteps参数，而非强制固定dt。

3. 对变量进行无量纲化，消除量级差异 #

如果原ODE中不同变量、系数的量级相差悬殊（比如有的是1e6，有的是1e-6），会导致数值计算精度丢失或溢出。此时建议无量纲化所有变量，把时间、状态变量都除以特征尺度M0，让所有量都处于O(1)量级：

# 无量纲化后：t' = t/M0，x' = x/M0
def dimensionless_ode(t_prime, x_prime):
    # 根据原ODE形式转换逻辑，示例如下：
    # 原方程 dx/dt = f(t, x) → 转换为 dx'/dt' = f(t'*M0, x'*M0)
    return your_transformed_expression(t_prime, x_prime)

无量纲化后再求解，能大幅降低数值量级带来的问题。

4. 调试原ODE函数的数值稳定性 #

求解失败也可能是原ODE函数内部计算出现异常（比如除以极小值、取对数负数等），导致出现NaN或inf，进而终止求解。建议在原ODE函数中加入调试输出：

def your_original_ode(t, x):
    # 打印中间值，排查是否出现异常数值
    print(f"当前t={t:.4f}, x={x}")
    # 原计算逻辑
    return result

这样可以定位到求解失败时的具体时间点和状态变量，快速找到问题根源。

总结：你的测试代码验证了当前量级下的步长设置是有效的，原问题的求解失败更可能是原ODE的刚性、求解器选型不当、变量量级差异过大，或者ODE函数内部的数值问题导致的，可以按照上述方向逐一排查。

内容的提问来源于stack exchange，提问作者Joshhh