I am trying to solve a set of coupled ODE, basically rate equations, all of them are of the form dy/dt = f(y), using python's scipy's odeint library. I am shocked by the time it is taking to solve it. I have solved the same set of equations before in FORTRAN, where each time step took about 0.1 sec. Here, the same time-step is taking 30 minutes! I know python is slower than FORTRAN or C++ but this seems impossible. What I have done: I have defined a function which returns a list, where each element is the value of dy_i/dt. Then I call odeint over a run it over small-time interval. The code is like this, if that helps.

def matrix_reac(abundances, time):  
    matrix = []
    densityfactor = gas_density(time)/gas_density(time0)
    for ii in range(len(abundances)):
        mat1, mat2 = 0.0, 0.0
        for ik in range(len(reac_rate)):
            reaction_rate = reac_rate[ik]*(gas_temperature(time)/300)**reac_nu[ik]*np.exp(-reac_actT[ik]/gas_temperature(time))
            if (species_list[ii] == network['R1'][ik] or species_list[ii] == network['R2'][ik]):
                if ('M' in [network['R1'][ik],network['R2'][ik]] or 'XX' in [network['R1'][ik],network['R2'][ik]]):
                    mat1 = mat1 - densityfactor*reaction_rate*abundances[ii]
                else: 
                    mat1 = mat1 - densityfactor**2*reaction_rate*abundances[np.where(network['R1'][ik]==species_list)]*abundances[np.where(network['R2'][ik]==species_list)]    
                    
            if (species_list[ii] == network['R4'][ik] or species_list[ii] == network['R5'][ik]):    
                if ('M' in [network['R1'][ik],network['R2'][ik]] or 'XX' in [network['R1'][ik],network['R2'][ik]]): 
                    mat2 = mat2 + densityfactor*reaction_rate*abundances[np.where(network['R1'][ik]==species_list)]
                else:
                    mat2 = mat2 + densityfactor**2*reaction_rate*abundances[np.where(network['R1'][ik]==species_list)]*abundances[np.where(network['R2'][ik]==species_list)]       
        
#        print (species_list[ii],mat1,mat2)
        matrix.append((mat1+mat2)*sec_1day)
    return matrix  

tstep = np.linspace(t_days[0],t_days[0]+1,2)
sol22 = odeint(matrix_reac,species_abun,tstep, atol=1.0e-6, rtol=1.0e-6)

I tried using mxstep, and it considerably shortened the time, but apparently the results are wrong. That is because if I increase the mxstep the yield is increasing rapidly with it, so I cannot make it small.
Could someone suggest something. In principle, I have to run it for a long time range, so if 1 time-step takes 30 mins, it will just be unfeasible. In my mind it seems something I am not doing right, it cannot be this slow, and 80 coupled ODE is not an absurdly large set anyway. Many thanks :)