Plotting bifurcation diagrams of Lorenz system with varying fractional order

Question

I need a Python code on how to plot bifurcation diagrams of Lorenz system with varying fractional order.

The one I've been able to come up with is this image and I am not sure it is correct

The genOrd(start, end, step) helps to get an array of numbers from start to end with increasing with step amount

The method for solving Fractional order Lorenz used here is the fractional order Adams-Bashforth-Moulton(ABM) numerical method for Fractional differential equations

x0 is initial condition = [-2.5, 6.6, -15.0]

def x_fn(x, y, z):
  return sigma * (y - x)
def y_fn(x, y, z):
  return (x * (rho - z)) - y
def z_fn(x, y, z):
  return (x*y) - (beta * z)

def Gamma(x):
  return gamma(x)


def ABM(alpha, N=10000, h=0.004, s=x0):
  b = np.zeros(N + 1)
  a = np.zeros(N + 1)
  x = np.zeros(N + 1)
  y = np.zeros(N + 1)
  z = np.zeros(N + 1)

  x[0] = s[0]
  y[0] = s[1]
  z[0] = s[2]

  for k in range(1, N + 1):
    a[k] = ((k+1)**(alpha + 1)) - (2 * (k**(alpha + 1))) + ((k-1)**(alpha + 1))
    b[k] = (k**alpha) - ((k-1)**alpha)

  def P(j, index, fn):
    factor = (h**alpha) / Gamma(alpha + 1)
    summ = 0
    for k in range(j+1):
      summ += b[j - k] * fn(x[k], y[k], z[k])
    if index == 0:
      return x[0] + (factor * summ)
    elif index == 1:
      return y[0] + (factor * summ)
    else:
      return z[0] + (factor * summ)

  def A(j, fn):
    summ = 0
    for k in range(1, j+1):
      summ += a[j - k] * fn(x[k], y[k], z[k])
    return summ

  def Not(j, fn):
    return (((j-1)**(alpha + 1)) - ((j - 1 - alpha) * (j**alpha))) * fn(x[0], y[0], z[0])

  F = (h**alpha) / Gamma(alpha + 2)

  for j in range(1, N+1):
    x_p = P(j, 0, x_fn)
    y_p = P(j, 1, y_fn)
    z_p = P(j, 2, z_fn)

    itx = Not(j, x_fn)
    ity = Not(j, y_fn)
    itz = Not(j, z_fn)

    x[j] = x[0] + F * (x_fn(x_p, y_p, z_p) + itx + A(j, x_fn))
    y[j] = y[0] + F * (y_fn(x_p, y_p, z_p) + ity + A(j, y_fn))
    z[j] = z[0] + F * (z_fn(x_p, y_p, z_p) + itz + A(j, z_fn))

  return x, y, z

The bifurcation codes

def getAndPlotBif():
  Z = []
  ordss = genOrds(0.8, 1.0, 0.001)
  ordz = []
  xs0, ys0, zs0 = x0
  for ord in ordss:
    x, y, z = ABM(ord, 3000, 0.004)
    for i in range(len(z)):
      Z.append(x[i])
      ordz.append(ord)
    xs0, ys0, zs0 = x[-1], y[-1], z[-1]
  Z = np.array(Z)
  ordz = np.array(ordz)
  return ordz, Z

This is the plotter code here:

def bif_plot(ords, Z, title="Bifurcation diagram of Lorenz system with different derivative order", x='Derivative order'):

  pl = plt

  pl.figure(figsize=(10, 6))

  pl.scatter(ords, Z, color='blue', s=0.5, alpha=0.2)

  pl.title(title)

  pl.xlabel(x)

  pl.ylabel("Z values")

  pl.show()

This is the result image