Plot basin of attraction for unforced duffing equation with matplotlib

I am evolving the system using a symplectic integrator of sixth order, this system has a pair of attractors on each side of the phase space. My attempt to plot the basins of attraction is by defining a grid of points in the phase space and evolving each of the points using the integrator, when the system enters a specific zone of the phase space (the attractor)(and what I can't manage to write correctly in the code) the integrator will stop and give a value 0 or 1 to the initial point, this will result in a grid of points with values 0 or 1, and when I plot the function's contour this will hopefully show the basins of attraction:

#yoshida 6th order symplectic integrator for attractor
import math as m
import numpy as np
import matplotlib.pylab as plt

#Hamilton equations
    def A(m0, a0, b0, x0, v0, gamma0): #p point
    return a0*x0/m0 - b0*x0*x0*x0/m0 -gamma0*v0
def V(m0,p0): #x point
    return p0/m0
def f(x0, p0, m0, a0, b0): #contour function
    return p0**2/(2*m0)-a0*x0**2/2+b0*x0**4/4


gamma=0.1 #dissipation parameter

#constants for integrator
dos=2.**(1./3.)
m=1.
a=0.25
b=0.01
a1=0.78451361047756 
b1=0.39225680523878 
a2= 0.23557321335936
b2=0.51004341191846
a3=-1.1776799841789
b3=-0.47105338540976 
a4=1.3151863206839
b4=0.068753168252520 
a5=-1.1776799841789
b5=0.068753168252520
a6=0.23557321335936
b6=-0.47105338540976
a7=0.78451361047756
b7=0.51004341191846
a8=0
b8=0.39225680523878

xc = np.linspace(-9, 9, 100) #contour grid
yc = np.linspace(-4, 4, 100) #contour grid

xp = np.linspace(-9, 9, 100) #grid for points
yp = np.linspace(-4, 4, 100) #grid for points   



XC, YC = np.meshgrid(xc, yc) #grid
ZC = f(XC, YC,m,a,b) 


def yoshi(xx,pp):  #function of yoshida integrator

    x=[xx] #initial position
    p=[pp] #initial momentum


    dt=0.001 #paso del integrador
    t=[0.] #lista de tiempo
    H=[p[0]*p[0]/(2.*m)-a*x[0]*x[0]/2.+b*x[0]*x[0]*x[0]*x[0]/4.+gamma*p[0]*x[0]] #list of the energy



    i=0
    while True: 
        e=2
        print (i)

        vini=p[i]/m

        v11=V(m,p[i])+b1*A(m,a,b,x[i],vini,gamma)*dt


        x11=x[i]+a1*v11*dt 


        v12=v11+b2*A(m,a,b,x11,v11,gamma)*dt


        x12=x11+a2*v12*dt


        v13=v12+b3*A(m,a,b,x12,v12,gamma)*dt


        x13=x12+a3*v13*dt


        v14=v13+b4*A(m,a,b,x13,v13,gamma)*dt

        x14=x13+a4*v14*dt


        v15=v14+b5*A(m,a,b,x14,v14,gamma)*dt

        x15=x14+a5*v15*dt


        v16=v15+b6*A(m,a,b,x15,v15,gamma)*dt


        x16=x15+a6*v16*dt


        v17=v16+b7*A(m,a,b,x16,v16,gamma)*dt


        x17=x16+a7*v17*dt


        v18=v17+b8*A(m,a,b,x17,v17,gamma)*dt


        x18=x17+a8*v18*dt


        x.append(x18)


        p.append(m*v18) 

        xi=x[i]
        pi=p[i]

        t.append(t[i]+dt)
        H.append(pi*pi/(2.*m)-a*xi*xi/2.+b*xi*xi*xi*xi/4.+gamma*pi*xi)
        if (pi<'0.86' and pi>'-0.86' and xi<'-3.5' and xi>'-6.17'): #my attempt for giving values
        e=0
            break
        if (pi<'0.86' and pi>'-0.86' and xi>'3.5' and xi<'6.17'):
            e=1
            break
        i=i+1
    if (e=='0'):
        return 0
    if (e=='1'):
        return 1



XP, YP = np.meshgrid(xp, yp) #grid for contour
ZP = yoshi(XP, YP) #function applied to the whole grid

#Attractors plot
plt.contourf(XP, YP, ZP, 2) 
plt.colorbar();
plt.title('Atractores')
plt.savefig('atractor.png', bbox_inches='tight')
plt.clf()

This attempt of the code seems to be trapped in the while loop forever, and I think it does not read correctly the "if" condition I am trying to impose. Is there a way I can solve this?