使用RK4方法求解非线性微分方程组时Decimal库引发decimal.InvalidOperation错误的排查求助

各位好，我现在正在用RK4方法求解一组非线性微分方程组，其中有个参数的量级达到了10^-121，为了保证计算精度，我引入了Decimal库，但现在一直遇到decimal.InvalidOperation的报错，折腾了好久都找不到问题根源，实在有点崩溃，恳请大家帮忙看看！

我的代码如下：

import numpy as np
from decimal import Decimal

def V(x,y,N):
    v = np.divide(y**2, (1-x**2-y**2)) * 3* H02 * (Om1 * np.exp(-3*N) + Om2 * np.exp(-4*N))
    return v

def f1 (x,y,l,q):
    return -3*x + l * np.sqrt(3/2) * y**2+ 3/2 * x * ( 2*x**2 + q*(1 - x**2 y**2))

def f2 (x,y,l,q):
    return - l * np.sqrt(3/2) * y * x + 3/2 * y * (2 * x**2 + q*(1 - x**2 - y**2))

def f3 (N, x, y, l, g, M, al):
    exp1,   G,   m  = Decimal(-2*g).exp(),    Decimal(g),         Decimal(M)
    xx,  a,   ll    = Decimal(x), Decimal(al),    Decimal(l)
    VV              = Decimal(V(x,y,N))
    log             = Decimal(VV/(m*exp1) + 1).ln()
    Tanh            = 1/G * log  - 1
    # gamm            = -(G*(1 - Tanh**2) - 2 * Tanh)/(np.sqrt(6*a)*ll)
    return  float(ll *xx/np.sqrt(a)* (G*(1 - Tanh**2) - 2 * Tanh + np.sqrt(6*a)*ll))

def R_K_4(xi,yi,li, g, M, al):
    N = int(round((lnaf-lnai)/dlna))
    lna = np.linspace(lnai, lnaf, N)
    x = np.zeros(N)
    y = np.zeros(N)
    l = np.zeros(N)
    x[0],y[0],l[0] = xi,yi,li
    for i in range(0,N-1):
        q    = 1 + 1/(3*(1 + 5780*np.exp(lna[i])))
        kx1 =  f1( x[i], y[i], l[i], q )
        ky1 =  f2( x[i], y[i], l[i], q )
        kl1 =  f3( lna[i], x[i], y[i], l[i], g, M, al)
        kx2 =  f1(                     x[i]+kx1*dlna/2,  y[i]+ky1*dlna/2,  l[i]+kl1*dlna/2, q )
        ky2 =  f2(                     x[i]+kx1*dlna/2,  y[i]+ky1*dlna/2,  l[i]+kl1*dlna/2, q )
        kl2 =  f3( lna[i] +  dlna/2,   x[i]+kx1*dlna/2,  y[i]+ky1*dlna/2,  l[i]+kl1*dlna/2,       g, M, al )
        kx3 =  f1(                     x[i]+kx2*dlna/2,  y[i]+ky2*dlna/2,  l[i]+kl2*dlna/2, q )
        ky3 =  f2(                     x[i]+kx2*dlna/2,  y[i]+ky2*dlna/2,  l[i]+kl2*dlna/2, q )
        kl3 =  f3( lna[i] +  dlna/2,   x[i]+kx2*dlna/2,  y[i]+ky2*dlna/2,  l[i]+kl2*dlna/2,        g, M, al )
        kx4 =  f1(                     x[i]+kx3*dlna,    y[i]+ky3*dlna,    l[i]+kl3*dlna,   q )
        ky4 =  f2(                     x[i]+kx3*dlna,    y[i]+ky3*dlna,    l[i]+kl3*dlna,   q )
        kl4 =  f3( lna[i] +  dlna,     x[i]+kx3*dlna,    y[i]+ky3*dlna,    l[i]+kl3*dlna,          g, M, al )
        x[i+1]  = x[i] +dlna/6*(kx1 + 2*kx2 + 2*kx3 + kx4)
        y[i+1]  = y[i] +dlna/6*(ky1 + 2*ky2 + 2*ky3 + ky4)
        l[i+1]  = l[i] +dlna/6*(kl1 + 2*kl2 + 2*kl3 + kl4)
    return x,y,l

当我用以下参数调用求解器时：

H02 = 1.44*1e-122
Om1, Om2 = 0.31, 5.45*1e-5
lnai,  lnaf, dlna  = -15, 10, 1e-2
x11, y11, l11 = R_K_4(0, 2.4e-11, -0.919, 128, 1e-10, 7/3)

就会触发decimal.InvalidOperation错误。我猜测错误可能在更早的计算步骤就已经产生了，但一直找不到具体位置，希望各位能帮我分析一下问题出在哪，谢谢大家！

备注：内容来源于stack exchange，提问作者user280016