title: Generating and Solving Simultaneous ODE's in Python
body: <p>I'm relatively new to Python, and am encountering some issues in writing a piece of code that generates and then solves a system of differential equations. </p>
  
  <p>My approach to doing this was to create a set of variables and coefficients, <em>(x0, x1, ..., xn)</em> and <em>(c0, c1 ,..., cn)</em> repsectively, in a list with the function <strong>var()</strong>. Then the equations are constructed in <strong>EOM1()</strong>. The initial conditions, along with the set of equations, are all put together in <strong>EOM2()</strong> and solved using <code>odeint</code>. </p>
  
  <p>Currently the code below runs, albeit not efficiently the reason for which I believe is because <code>odeint</code> runs through all the code with each interaction (that's something else I need to fix but isn't the main problem!).</p>
  
  <pre><code>import sympy as sy
  from scipy.integrate import odeint
  
  n=2
  cn0list = [0.01, 0.05]
  xn0list = [0.01, 0.01]
  
  def var():
      xnlist=[]
      cnlist=[]
      for i in range(n+1):
          xnlist.append('x{0}'.format(i))
          cnlist.append('c{0}'.format(i))
      return xnlist, cnlist
  
  def EOM1():
      drdtlist=[]
      for i in range(n):
          cn1=sy.Symbol(var()[1][i])
          xn0=sy.Symbol(var()[0][i])
          xn1=sy.Symbol(var()[0][i+1])
          eom=cn1*xn0*(1.0-xn1)-cn1*xn1-xn1
          drdtlist.append(eom)
      xi=sy.Symbol(var()[0][0])
      xf=sy.Symbol(var()[0][n])
      drdtlist[n-1]=drdtlist[n-1].subs(xf,xi)
      return drdtlist
  
  def EOM2(xn, t, cn):
      x0, x1 = xn
      c0, c1 = cn
      f = EOM1()
      output = []
      for part in f:
          output.append(part.evalf(subs={'x0':x0, 'x1':x1, 'c0':c0, 'c1':c1}))
      return output
  
  abserr = 1.0e-6
  relerr = 1.0e-4
  stoptime = 10.0
  numpoints = 20
  
  t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]
  
  wsol = odeint(EOM2, xn0list, t, args=(cn0list,), atol=abserr, rtol=relerr)
  </code></pre>
  
  <p>My problem is that I had difficulty getting Python to treat the variables generated by Sympy appropriately. I got around this with the line</p>
  
  <pre><code>output.append(part.evalf(subs={'x0':x0, 'x1':x1, 'c0':c0, 'c1':c1}))
  </code></pre>
  
  <p>in <strong>EOM2()</strong>. Unfortunately, I do not know how to generalize this line to a list of variables from <em>x0</em> to <em>xn</em>, and from <em>c0</em> to <em>cn</em>. The same applies to the earlier line in <strong>EOM2()</strong>,</p>
  
  <pre><code>    x0, x1 = xn
      c0, c1 = cn
  </code></pre>
  
  <p>In other words I set <em>n</em> to an arbitrary number, is there a way for Python to interpret each element as it does with the ones I manually entered above? I have tried the following</p>
  
  <pre><code>output.append(part.evalf(subs={'x{0}'.format(j):var(n)[0][j], 'c{0}'.format(j):var(n)[1][j]}))
  </code></pre>
  
  <p>yet this yields the error that led me to use <strong>evalf</strong> in the first place,</p>
  
  <pre><code>TypeError: can't convert expression to float
  </code></pre>
  
  <p>Is there any way do what I want to, generate a set of <em>n</em> equations which are then solved with <code>odeint</code>?</p>