Python simulation of the Van der Pol Oscillator

One of the things I'm doing at the moment is watching an excellent series of MIT OCW lectures by Assoc. Prof. Russ Tedrake on underactuated robotics. I've learned a lot watching the videos, which pick up at a point where some parts of my doctoral studies finished.

Just to pick something small that is at least tangentially related to the related to the course, I've attached a python snippet of code for generating the Van der Pol Oscillator. This system is significant in that it exhibits limit cycles, which are an important tool for reasoning about walking robots

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

mu = 0.2

def van_der_pol_oscillator_deriv(x, t):
    nx0 = x[1]
    nx1 = -mu * (x[0] ** 2.0 - 1.0) * x[1] - x[0]
    res = np.array([nx0, nx1])
    return res

ts = np.linspace(0.0, 50.0, 500)

xs = odeint(van_der_pol_oscillator_deriv, [0.2, 0.2], ts)
plt.plot(xs[:,0], xs[:,1])
xs = odeint(van_der_pol_oscillator_deriv, [-3.0, -3.0], ts)
plt.plot(xs[:,0], xs[:,1])
xs = odeint(van_der_pol_oscillator_deriv, [4.0, 4.0], ts)
plt.plot(xs[:,0], xs[:,1])
plt.gca().set_aspect('equal')
plt.savefig('vanderpol_oscillator.png')
plt.show() 

Which when simulated looks something like this: