def METHOD_NAME(system, f):
    def _func(t, rho, dt=1e-6):
        n = rho.shape[0]
        f0 = f(t, rho)
        out = (f(t+dt, rho) - f0) / dt
        jac = np.zeros((n, n), dtype=complex)
        for i in range(n):
            dxi = basis(n, i).data
            jac[:, i] = (f(t, rho + dt * dxi) - f0).to_array().flatten() / dt
        out = out + _data.Dense(jac) @ system.drift(t, rho)
        for i, j in product(range(n), repeat=2):
            dxi = basis(n, i).data
            dxj = basis(n, j).data
            sec = f(t, (rho + dxi * dt + dxj * dt))
            sec = sec - f(t, (rho + dxj * dt))
            sec = sec - f(t, (rho + dxi * dt))
            sec = sec + f0
            sec = sec / dt / dt * 0.5
            for k in range(system.num_collapse):
                out = out + (
                    sec
                    * _data.inner(dxi, system.diffusion(t, rho)[k])
                    * _data.inner(dxj, system.diffusion(t, rho)[k])
                )
        return out
    return _func