So, I've seen the coded solution to my question in Mathematica, but with very little understanding of mathematica, I havn't been able to reproduce it yet.
This is what I'm trying to do with Python: https://mathematica.stackexchange.com/questions/159211/how-to-make-a-bifurcation-diagram-of-the-lorenz-system-under-a-varying-parameter
I'm thinking my errors are in understanding how to calculate what I'm looking for and how to adjust my visualization to make it look just like that in the link, but any ideas are welcome.
The code I have so far looks like this:
def lorenz_system(x,y,z,r,s=10,b=6):
x_dot = s*(y-x)
y_dot = r*x-y-x*z
z_dot = x*z-b*z
return x_dot, y_dot, z_dot
dr = 0.1 # parameter step size
r=np.arange(40,200,dr) # parameter range
dt = 0.001 # time step
t = np.arange(0,10,dt) # time range
#initialize solution arrays
xs = np.empty(len(t) + 1)
ys = np.empty(len(t) + 1)
zs = np.empty(len(t) + 1)
#initial values x0,y0,z0 for the system
xs[0], ys[0], zs[0] = (1, 1, 1)
for R in r:
for i in range(len(t)):
#approximate numerical solutions to system
x_dot, y_dot, z_dot = lorenz_system(xs[i], ys[i], zs[i],R)
xs[i + 1] = xs[i] + (x_dot * dt)
ys[i + 1] = ys[i] + (y_dot * dt)
zs[i + 1] = zs[i] + (z_dot * dt)
#calculate and plot the peak values of the z solution
for i in range(0,len(zs)-1):
#using only the positive values in the z solution
if zs[i]>0:
#find the local maxima
if (zs[i-1] < zs[i] and zs[i] > zs[i+1]):
if (zs[i]<=1000):
#plot the local maxima point of the z solution that used the parameter R in r
plt.scatter(R,zs[i], color='black')
plt.xlim(0,200)
plt.ylim(0,400)