solve_ivp求解phi/psi方程组耗时久、结果过大问题求助

我正在使用以下代码求解phi、psi对应的方程组：

def _derivs(self, t, y):
    """
    定义待求解微分方程的函数：
     1. I * phi_ddot = tau_motor - tau_drag.
     2. I * psi_ddot = F_tot x R - kappa * psi - gamma * psi_dot
    """

    phi, phi_dot, psi, psi_dot = y
    phi_ddot = 0.02 - phi_dot * np.abs(phi_dot)
    f_lift = 1.7 * phi_dot * np.abs(phi_dot)
    f_drag = - 0.8 * phi_dot * np.abs(phi_dot)
    f_tot = np.sqrt(f_drag ** 2 + f_lift ** 2)
    psi_ddot = f_tot - psi - psi_dot
    dy_dt = [phi_dot, phi_ddot, psi_dot, psi_ddot]
    return dy_dt

phi0, phi_dot0, psi0, psi_dot0= 0, 2e-4, 0, 2e-4
start_t, end_t = 0,0.2
sol = solve_ivp(_derivs,t_span=(start_t, end_t),y0=[phi0, phi_dot0, psi0, psi_dot0])

注：代码中的1.7、-0.8等常数前缀无需关注（实际计算方式不同，与当前问题无关）

当前遇到的问题：求解器耗时数分钟才完成，返回的sol.y结果长度不合理（超过10万条数据）；但仅求解phi时，求解器返回的结果长度正常。

• 系统强非线性引发刚度问题：psi的微分方程依赖f_tot，而f_tot是phi_dot的平方项函数，属于强非线性项。当phi_dot变化时，psi_ddot会出现突变式的非线性变化，这会让求解器判定系统存在刚度特性。为了满足误差控制要求，求解器会自动把积分步长缩到极小，直接导致计算量暴增、结果数据量异常庞大。单独求解phi时，方程的非线性更平缓，求解器能维持正常步长。
• 默认求解器不适配强非线性系统：solve_ivp默认采用RK45求解器，这类显式求解器对强非线性、带刚度的系统适应性很差。引入psi的强非线性方程后，RK45需要不断细分步长来跟踪解的快速变化，导致计算效率急剧下降。而单独求解phi的方程非线性程度低，RK45可以高效处理。
• 初始微小值的非线性放大：初始条件中phi_dot0和psi_dot0都是2e-4量级的极小值，经过平方项计算后，f_tot的变化会被非线性放大，使得psi的方程出现快速波动的解。求解器为了精准捕捉这种快速变化，不得不持续使用极小步长积分，最终导致数据量异常。

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