Bifurcation diagram

code

import math
import matplotlib.pyplot as plt

def Df(LIM):
    w=[0 for x in range(0,LIM)]
    theta=[0 for x in range(0,LIM)]
    t=[0 for x in range(0,LIM)]
    w[0]=0;theta[0]=0.2;Omega_D = 2/3;dt=0.01;g=9.8;l=9.8;q=0.5
    P_theta=[]
    P_FD=[]
    for j in range(0,150):
        F_D=1.35+0.001*j

        for i in range(0,LIM-1):
            w[i+1]=w[i] + (-(g/l) * math.sin(theta[i]) - q*w[i] + F_D * math.sin(Omega_D * t[i]) )*dt
            theta[i+1]=theta[i]+w[i+1]*dt
            if theta[i+1] > math.pi:
                theta[i+1]=theta[i+1]-2*math.pi
            elif theta[i+1] < -math.pi:
                theta[i+1]=theta[i+1]+2*math.pi
            else:
                theta[i+1]=theta[i+1]
            t[i+1]=t[i]+dt
            if (t[i+1]>900*math.pi)and((Omega_D*t[i+1])/(2*math.pi)-int((Omega_D*t[i+1])/(2*math.pi))<(Omega_D*dt)/(2*math.pi)):
                P_theta.append(theta[i+1])
                P_FD.append(F_D)
    return [ P_FD , P_theta]
plt.scatter(Df(350000)[0],Df(350000)[1],s=1)