title: Vectorize Forward Euler method for system of differential equations
body: <p>I am numerically solving for x(t) for a system of first order differential equations. The system is:</p>
  
  <p>dx/dt = y  <br />
  dy/dt = -x - a*y(x^2 + y^2 -1)</p>
  
  <p>I have implemented the Forward Euler method to solve this problem as follows:</p>
  
  <pre><code>def forward_euler():
      h = 0.01
      num_steps = 10000
  
      x = np.zeros([num_steps + 1, 2]) # steps, number of solutions
      y = np.zeros([num_steps + 1, 2])
      a = 1.
  
      x[0, 0] = 10. # initial condition 1st solution
      y[0, 0] = 5.
  
      x[0, 1] = 0.  # initial condition 2nd solution
      y[0, 1] = 0.0000000001
  
      for step in xrange(num_steps):
          x[step + 1] = x[step] + h * y[step]
          y[step + 1] = y[step] + h * (-x[step] - a * y[step] * (x[step] ** 2 + y[step] ** 2 - 1))
  
      return x, y
  </code></pre>
  
  <p>Now I would like to vectorize the code further and keep x and y in the same array, I have come up with the following solution:</p>
  
  <pre><code>def forward_euler_vector():
      num_steps = 10000
      h = 0.01
  
      x = np.zeros([num_steps + 1, 2, 2]) # steps, variables, number of solutions
      a = 1.
  
      x[0, 0, 0] = 10. # initial conditions 1st solution
      x[0, 1, 0] = 5.  
  
      x[0, 0, 1] = 0.  # initial conditions 2nd solution
      x[0, 1, 1] = 0.0000000001
  
      def f(x): 
          return np.array([x[1],
                           -x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1)])
  
      for step in xrange(num_steps):
          x[step + 1] = x[step] + h * f(x[step])
  
      return x
  </code></pre>
  
  <p>The question: forward_euler_vector() works, but was this to best way to vectorize it? I am asking because the vectorized version runs about 20 ms slower on my laptop:</p>
  
  <pre><code>In [27]: %timeit forward_euler()
  1 loops, best of 3: 301 ms per loop
  
  In [65]: %timeit forward_euler_vector()
  1 loops, best of 3: 320 ms per loop
  </code></pre>