Avoiding plotting ODEs divergent solutions ODEint

Question

I'm trying to plot a phase plane and I want it to look nice. However, some solutions of the system of equations diverge because of the initial conditions. Is there some way that I can make a try/except chain in order when the solution diverges it doesn't plot it. Here is my code:

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

def aux_func(x):  
  y = x[0]-x[1]
  if (np.abs(y) <= 1):
    f = y**3 + 0.5*y 
  else: 
    f = 2*y - np.sign(y)
  return f 

def function(x,t): 
  x1_dot = x[1]  
  x2_dot = -x[1] - aux_func(x)  
  return [x1_dot,x2_dot]

ts = np.linspace(0, 20, 300)  
ic_1 = np.linspace(-1,1,10)
ic_2 = np.linspace(-1,1,10) 

for r1 in ic_1:
  for r2 in ic_2:
    x0 = (r1,r2)
    try:
      xs = odeint(function, x0, ts)
      plt.plot(xs[:,0], xs[:,1],"r-",linewidth=.8) 
    except: 
      pass

# Nombre de los ejes, limites,
plt.xlabel("$x_1$", fontsize=12)
plt.ylabel("$x_2$", fontsize=12)
# plt.tick_params(labelsize=10)
# plt.xticks(np.linspace(0,1,11))
# plt.yticks(np.linspace(0,1,11))
plt.xlim(-1, 1)
plt.ylim(-1, 1)
# Grafica el campo vectorial
X1, X2 = np.mgrid[-1:1:20j,-1:1:20j]
u=X2
d= X1-X2

t = np.zeros(np.shape(d))
for i in range(len(d)):
    for j in range(len(d[0])):
        if np.abs(d[i][j]) > 1:
            t[i][j]= 2*d[i][j]-0.5*np.sign(d[i][j])
        else:
            t[i][j] =d[i][j]**3 + 0.5*d[i][j]
v=-X2-t

pl.quiver(X1, X2, u, v, color = 'b',width = .002)
plt.grid()
plt.title('Plano de Fase Punto 1')
#plt.savefig('FasePunto4.png')
plt.show()

The code is plotting the following:

Appreciate the help.

Answer 1

This can be solved by avoiding the wrong divergences at all, so that there is no need for exception handling.

This is a discontinuous ODE which can lead to unusual effects like a sliding mode. One way to quickly work around that is to mollify the jump by implementing a blending zone where the vector field changes quickly but continuously from one phase to the other (see Unsure about how to use event function in Matlab for other generic work-arounds). The changes for that can be implemented as

def aux_func(x):  
  def softsign(u): return np.tanh(1e4*u)
  y = x[0]-x[1]
  h = 0.5*(1+softsign(y**2-1)
  # h is about zero for |y|<1 and about 1 for |y|>1
  f1 = y**3 + 0.5*y          # for |y|<1
  f2 = 2*y - softsign(y)     # for |y|>1, note the second mollification
  return (1-h)*f1+h*f2 

With no further changes to the code this gives the plot

Note that pylab is obsolete, all its functionality can also be accessed via plt=matplotlib.pyplot.