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?