title: How to keep a consistent order to eigenvalues of a matrix using numpy or scipy?
body: <p>I am trying to get a plot of a set of eigenvalues of a matrix against a detuning factor epsilon. I want the plot to look like this
  <a href="https://i.stack.imgur.com/MPVgT.png" rel="nofollow noreferrer">Matlab Plot using eig()</a></p>
  <p>However when I use <code>np.linalg.eigvals</code> I get the following
  <a href="https://i.stack.imgur.com/mux6j.png" rel="nofollow noreferrer">Python eigvals plot</a></p>
  <p>I also tried using <code>np.linalg.eigvalsh</code> which gave <a href="https://i.stack.imgur.com/jo4Ws.png" rel="nofollow noreferrer">Python eigvalsh plot</a>.</p>
  <p>The problem seems to be how the eigenvalues are ordered upon the return for the function. I was wondering if there's any way to get it so I produce plot lines like in the first image from matlab. I've also tried the equivalent scipy functions which just gave the same as the numpy.</p>
  <p>Here is a copy of my code</p>
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  
  mu = 57.88e-3
  eps = np.linspace(-0.26,0.26,100)
  
  def Model_g1g2(g1, g2, B, t):
     
      E = np.empty([len(eps), 5])  # Initialise array for eigenvalues
      for i in range(len(eps)):
          # Matrix A
          A = np.array([[(-eps[i]+(g1+g2)*mu*B)/2, 0, 0, 0, 0], 
                    [0, -eps[i]/2, 0, (g1-g2)*mu*B, 0],
                    [0, 0, (-eps[i]-(g1+g2)*mu*B)/2, 0, 0],
                    [0, (g1-g2)*mu*B/2, 0, -eps[i]/2, t],
                    [0, 0, 0, t, eps[i]/2]
                    ])
  
          E[i,:] = np.linalg.eigvals(A) # Store eigenvalues
      return E
  
  E = Model_g1g2(1, 4, 0.5, 0.06)
  
  # Produce Plot
  for i in range(5):
      plt.plot(eps, np.real(E[:,I]))
  
  plt.show()
  </code></pre>