I have data points to fit with a model. However, my model is not difined as an analytic function, but as the numerical solution of a ODE system

Question

Step 1: solving edo system

#Defining EDO's systems

def H(w, a, Omega_m0, H0, M, epsilon, b, g):   
    H, x, y = w 
            ...
    dwda = [x, y, "function of the x and y - f(x,y)"] 
    return dwda

#Initial conditions
w0 = [x0, y0, f(0)]

#Independent variable
a = np.linspace(0.2, 1, 1000)

#solution
sol = odeint(H, w0, a, args=(Omega_m0, H0, M, epsilon, b, g))

Step 2: Fitting numerical model to the data points

#Loading of data points
data = np.loadtxt('data.txt')
z = data[:, 0]
Hz = data[:, 1]
Erro_Hz = data[0:, 2]
plt.scatter(z, Hz)
plt.errorbar(z, Hz, Erro_Hz, linestyle='None')

#Defining numerical solution as the model

def Func(z, Omega_m0, H0, M, epsilon, b, g):
    return sol[:,0]

fmodel=Model(Func)
result = fmodel.fit(Hz, z=z, Omega_m0=??, H0=??, M=??, epsilon=??, b=??, g=??)
print(result.fit_report())

Step 3: I need to connect step 1 with step 2 and get the best Omega_m0, H0, M, epsilon, b, g.