title: Nullcline Plot for Nonlinear System of ODEs
body: <p>I am attempting to plot the nullcline (steady state) curves of the Oregonator model to assert the existence of a limit cycle by applying the Poincare-Bendixson Theorem. I am close, but for some reason the plot that is produced shows two straight lines. I think it has something to do with the plotting stage. Any ideas?</p>
  
  <p>Also any hints for how to construct a quadrilateral to apply the theorem with would be most appreciated. </p>
  
  <p>Code: </p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  
  # Dimensionless parameters
  eps = 0.04
  q = 0.0008
  f = 1
  
  
  # Oregonator model as numpy array
  def Sys(Y, t = 0):
      return np.array((Y[0] * (1 - Y[0] - ((Y[0] - q) * f * Y[1]) / (Y[0] + q)) / eps, Y[0] - Y[1] ))
  
  
  
  # Oregonator model steady states
  def g(x,z):
      return (x * (1 - x) + ((q - x) * f * z) / (q + x)) / eps
  def h(x,z):
      return x - z
  
  # Initial lists containing values
  x = []
  z = []
  
  def sys(iv1, iv2, dt, time):
      # initial values:
      x.append(iv1)
      z.append(iv2)
      # Compute and fill lists
      for i in range(time):
          x.append(x[i] + (g(x[i],z[i])) * dt)
          z.append(z[i] + (h(x[i],z[i])) * dt)
      return x, z
  
  sys(1, 0.5, 0.01, 30)
  
  # Locate and find equilibrium points
  eqp = []
  
  def find_fixed_points(r):
      for x in range(r):
          for z in range(r):
              if ((g(x, z) == 0) and (h(x, z) == 0)):
                  eqp.append((x,z))
  
      return eqp
  
  
  # Plot nullclines
  plt.plot([0,2],[2,0], 'r-', lw=2, label='x-nullcline')
  plt.plot([1,1],[0,2], 'b-', lw=2, label='z-nullcline')
  
  # Plot equilibrium points
  for point in eqp:
      plt.plot(point[0],point[1],"red", marker = "o", markersize = 10.0)
  
  plt.legend(loc='best')
  
  x = np.linspace(0, 2, 20)
  z = np.linspace(0, 2, 20)
  
  X1 , Z1  = np.meshgrid(x, z)                    # Create a grid
  DX1, DZ1 = Sys([X1, Z1])                        # Compute reaction rate on the grid
  M = (np.hypot(DX1, DZ1))                        # Norm reaction rate
  M[ M == 0] = 1.                                 # Avoid zero division errors
  DX1 /= M                                        # Normalise each arrows
  DZ1 /= M
  
  plt.quiver(X1, Z1, DX1, DZ1, M, pivot='mid')
  plt.xlabel("x(\u03C4)")
  plt.ylabel("z(\u03C4)")
  plt.legend()
  plt.grid()
  plt.show()
  </code></pre>