title: Odeint for solving 80 couple differential equations
body: <p>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.</p>
  <pre><code>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)
  </code></pre>
  <p>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.<br />
  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 :)</p>