title: Solving a set of ODE with SciPy
body: <p>I'm trying to solve a set of ODE to simulate starch hydrolysis by amylases (enzymes). When I try to solve the set of equation, I get a </p>
  
  <pre><code>lsoda--  at current t (=r1), mxstep (=i1) steps   
         taken on this call before reaching tout     
        In above message,  I1 =       500
        In above message,  R1 =  0.6333483931400E+00
  Excess work done on this call (perhaps wrong Dfun type).
  Run with full_output = 1 to get quantitative information.
  </code></pre>
  
  <p>error. The code :</p>
  
  <pre><code>#!/usr/bin/env python
  # -*- coding: utf-8 -*-
  
  import numpy as np
  import math
  from scipy.integrate import odeint
  
  def TemperatureProfile (timeandtemp,t) :
  
      if t==0.0 :
          T = 273.15+55
      elif t&gt;0.0 and t&lt;9.0 :
          T = 273.15+55+t
      elif t&gt;= 9.0 and t&lt;=9.0+timeandtemp[0][0] :
          T = 273.15+ timeandtemp[0][1]
      elif t&gt;9.0+timeandtemp[0][0] and t&lt;9.0+timeandtemp[0][0]+8.0 :
          T = 273.15+ timeandtemp[0][1]+t
      elif t&gt;=9.0+timeandtemp[0][0]+8.0 and t&lt;=9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0] :
          T = 273.15+timeandtemp[1][1]
      elif t&gt;9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0] and t&lt;9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0]+6:
          T = 273.15+ timeandtemp[1][1]+t
      elif t&gt;=9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0]+6 and  t&lt;9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0]+6 :
          T = 273.15+timeandtemp[2][1]
  
      return T
  
  def Flux (ci,temperature) :
  
      '''
          Constants
      '''
      K_gel1 = 5.7*10**31
      K_gel2 = 3.1*10**14
  
      k_glc = 0.023
      k_glc2 = 2.9*10**-8
      k_alphamal = 0.389
      k_betamal = 0.137
      k_alphamal2 = 1.2*10**-7
      k_betamal2 = 8.4*10**-8
      k_mlt = 0.117
      k_mlt2 = 1.5*10**-8
      k_dex = 0.317
      k_degalpha = 6.9*10**30 
      E_degalpha = 224.2
      k_degbeta = 7.6*10**60
      E_degbeta = 410.7
  
  
      starchUG = ci[0]
      starchG = ci[1]
      glc = ci[2]
      mal = ci[3]
      mlt = ci[4]
      dex = ci[5]
      enzA = ci[6]
      enzB = ci[7]
  
  
      '''
          Relative activities
      '''
      if temperature &lt; 336.15 :
          a_alpha = -0.0011*temperature**3 + 1.091*temperature**2 - 352.89*temperature + 38008.3
          a_beta = 0.049*temperature - 13.9 
      else :
          a_alpha = 0.0055*temperature**3 - 5.663*temperature**2+ 1941.9*temperature- 221864 
          a_beta = 0.374*temperature + 128.3 
  
      ''' 
          Equations
      '''
      if temperature &lt; 333.15 :
          r_gel = K_gel1*starchUG
      else :
          r_gel = K_gel2*starchUG
  
      r_glc = k_glc*a_alpha*glc
      r_mal = (k_alphamal*a_alpha+k_betamal*a_beta)*starchG
      r_mlt = k_mlt*a_alpha*starchG
      r_dex = k_dex*a_alpha*starchG
      r_glc2 = k_glc2*a_alpha*dex
      r_mal2 = k_alphamal2*a_alpha*dex+k_betamal2*a_beta*dex
      r_mlt2 = k_mlt2*a_alpha*dex
      r_degA = k_degalpha*math.exp(-E_degalpha/(8.3145*temperature))*enzA
      r_degB = k_degbeta*math.exp(-E_degbeta/(8.3145*temperature))*enzB
      r_acA = k_degalpha*math.exp(-E_degalpha/(8.3145*temperature))*a_alpha*enzA
      r_acB = k_degbeta*math.exp(-E_degbeta/(8.3145*temperature))*a_beta*enzB
  
      return np.array([r_glc,r_mal,r_mlt, r_dex, r_glc2, r_mal2, r_mlt2, r_degA, r_degB, r_acA, r_acB])
  
  
  
  
  def Secmembre (ci,t,tempProf) :
  
  
      temperature = TemperatureProfile(tempProf,t)
  
      mS = np.mat([
  [-1, 0, 0,  0,  0,  0,  0,  0,  0,  0,  0],
  [1, 0,  0,  0,  0,  0,  0,  0,  0,  0,  0], 
  [0, 1,  0,  0,  0,  1,  0,  0,  0,  0,  0], 
  [0, 0,  1,  0,  0,  0,  1,  0,  0,  0,  0], 
  [0, 0,  0,  1,  0,  0,  0,  1,  0,  0,  0], 
  [0, 0,  0,  0,  1,  -1, -1, -1, 0,  0,  0]])
  
      mS = np.transpose(mS)
      mS = np.array(mS)
      metabos = Flux(ci,temperature)
      varMetabos = np.dot(metabos,mS)
  
      return varMetabos
  
  
  
  # Initial conditions
  
  
  initCond = []
  
  temperatureProfile = ((30,64),(30,72),(5,78)) # Temperature steps in minutes
  
  t0 = 0
  tf  = temperatureProfile[0][0]+temperatureProfile[1][0]+temperatureProfile[2][0]+9+8+6
  nbPoints = 100
  
  timeProfile = np.linspace(t0,tf,nbPoints)
  print timeProfile
  initCond.append(113.5)
  initCond.append(0)
  initCond.append(4)
  initCond.append(5)
  initCond.append(0)
  initCond.append(0)
  initCond.append(80000)
  initCond.append(80000)
  
  initCond=np.array(initCond)
  
  result = odeint(Secmembre,initCond,timeProfile,args=(temperatureProfile,))
  </code></pre>
  
  <p>When working on Matlab I used to solve with ode23tb. I think odeint is not the right solver but I don't know which one to use. Maybe someone is familiar with this type of equation ?</p>