关于分岔图的一些问题，跪求解答

是一个非线性动力学方程求解并画出它的分岔图，在文献上画是这个样子的：

而我画出来的则是这个样子的虽然我只画了一部分，但是差别也太大了

这个是我画的1.0680到1.0690的分岔图，横坐标有1000个点，在文献上看应该是个二周期的，而我这个就很乱，希望有谁能教教我，救救孩子吧 😦

下面附上用到的代码

# coding=gbk
from scipy.integrate import odeint
import numpy as np
import math
import scipy.signal as signal
from sys import exit
#----------------------------------------------
#import time
#start = time.perf_counter()
#----------------------------------------------
import csv
import datetime
start=datetime.datetime.now()
print (start)

def math_model(w, t, D, N, B, r):
    # 给出位置矢量w，和三个参数v, k, r计算出
    # dx/dt, dy/dt, dz/dt的值
    x, y, z= w
    # 直接与math_model的计算公式对应
    return np.array([y, (-2)*B*y - np.power(N,2) * np.sin(x) + r * np.power(N,2) * np.cos(z), D])
    # K=k,m=w0,A=A,B=r,o=seita

t = np.arange(0, 50, 0.0001)  # 创建时间点#############################################################################


#调用ode对math_model进行求解, 用两个不同的初始值

#定义两种分岔图函数，给名字和r值，返回r值下的局部极值数组   D          N          B      r
def F_bifurcation(function_name, r):
    function_name = odeint(math_model, (1, 1, 1), t, args=( 2*np.pi , 3*np.pi , 3/4*np.pi, r))
    #输出的是一个三列，行数为t变换次数的矩阵，取第一列，第二列分别做分岔图，取列方法：a[:,1](第二列)。
    F1_column = function_name[:, 0]  # 取矩阵第一列
    F_column = F1_column[400000:500000]  #取第一列的后10000个数据点########################################################################
    # print (F_column[signal.argrelextrema(F_column, np.greater)])
    return F_column[signal.argrelextrema(F_column, np.greater)]  #返回后10000个数据点的相对极值


# 绘图
# from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

plt.figure(figsize=(20, 10))

e = np.arange(1.068, 1.069, 0.000001)   #对应r的取值数组##########################################################################################111
#  绘制矩阵第一列极值对应点的分岔图
#  两次循环，r确定以后得到对应极值点，作图，r增加一个步长以后得到对应极值点，作图，以此类推。


with open('test28.csv','a')as f:##############################################################################################################2222
    f_csv = csv.writer(f)
    for x in range(len(e)):
        print(x)
        d = F_bifurcation('nb' + str(x), e[x])
        addd = np.insert(d,0,e[x])
        RE = addd.reshape((1, len(d)+1))
        print(RE)
        f_csv.writerows(RE)
        #d = F_bifurcation('nb' + str(x), round(e[x], 5))
        for y in range(len(d)):
            plt.scatter(e[x], d[y], marker='.', edgecolor='none', linewidths=0.1, s=10)


plt.title("F_bifurcation", fontsize=24)
plt.xlabel("r", fontsize=14)
plt.ylabel("relative_extreme", fontsize=14)


#----------------------------------------
#end = time.perf_counter()
#print('Running time: %s Seconds' % (end - start))
#----------------------------------------
end=datetime.datetime.now()
print('Running time: %s Seconds'%(end-start))

x = np.arange(1.068, 1.069, 0.0001)# x横坐标##################################################################################################33333
plt.grid()
plt.xticks(x, rotation='vertical')
plt.xlim(1.068, 1.069)######################################################################################################################444444

plt.savefig('D:\python_work\Geany_work\damped_driven_pendulum_bif2\power/测试t500000_50_0.0001_r1 1.068-1.069-1000.eps')
plt.savefig('D:\python_work\Geany_work\damped_driven_pendulum_bif2\power/测试t500000_50_0.0001_r1 1.068-1.069-1000.svg')######################################
plt.savefig('D:\python_work\Geany_work\damped_driven_pendulum_bif2\power/测试t500000_50_0.0001_r1 1.068-1.069-1000.png')

plt.show()

# plt.xlim(0,10)

user_input_str = input('input string:')
if user_input_str == 'q':
    exit()
#45,50,81