Poincare Section + Bifurcation of Verlet Algorithm for chaotic oscillator

Question

I have been working on a verlet algorithm to give the solutions to damped driven oscillator problems. So far my code for the oscillator itself has worked and given sensible results for phase space and config space. I'm trying to produce poincare cuts and a bifurcation plot of the same code now, having trouble implementing it.

My code:

import numpy as np
import matplotlib.pyplot as plt


theta0 = np.pi/6. #start angle - vary this
w0 = 0. #start velocity
t0 = 0. #start time
gam = 0.5 #damping
Q = 1.5 #drive factor
nt = 500 #number of steps over time period
nst = 200*nt #number of periods
W = 2./3. #vary this

dt = (1/20*np.pi/W)/nt #time interval
T = int(nst*2*np.pi) #total time to plot for

#Start Arrays
t = np.array([t0])
theta = np.array([theta0])
w = np.array([w0])
theta_mid = np.array([])
w_mid = np.array([])
a_mid = np.array([])

#Defining Angular Accelation
def alp(theta,w,t):
    return  Q*np.sin(W*t)-gam*w-np.sin(theta)

#Potential energy
def pot(theta):
    return np.cos(theta)

#kinetic energy
def kin(w):
    return 0.5 * w**2

#Euler Richardson
for i in range(1,nst+1):
    theta_mid =np.append(theta_mid, theta[i-1]+w[i-1]*0.5*dt)
    w_mid =np.append(w_mid, w[i-1]+alp(theta[i-1],w[i-1],t[i-1])*0.5*dt) 
    a_mid =np.append(a_mid, alp(theta_mid[i-1], w_mid[i-1],t[i-1]+dt/2)) 
    theta=np.append(theta, theta[i-1]+w_mid[i-1] *dt)
    w=np.append(w,w[i-1]+a_mid[i-1]*dt)
    t = np.append(t,t[i-1]+dt)
    n = t[i]*W/2*np.pi
        
#The rest of the code is identical to part 2, we just plot space against velocity rather than time
plt.rc('font', family = 'serif', serif = 'cmr10',size=10) #to make like latex
plt.title(("Phase Space diagram for Q={}").format(Q))
plt.xlabel('Theta')
plt.ylabel('Omega')
#plt.plot(theta,w)
plt.scatter(theta[int(nst*0.5)::nt],w[int(nst*0.5)::nt])
plt.show()

Anyone got any ideas how to better implement it? Feeling like the time step and scatter part are not quite certain in my head. Some people have introduced an if loop in the main function but I dont see how it helps..

tried to plot for only multiples of 2pi plotted nothing or something that resembles the non stroboscopic