Bifurcation2D.sagews

#######################
#Bifurcation Diagram 2D
#by Sinan Kapcak
#2017
#www.k-interact.net/dds
#######################

a,x,y=var('a x y')
@interact
def _(
    f=input_box(label='$f(x,y)=$',default=a*y/(x^2 + y+1),width=60),
    g=input_box(label='$g(x,y)=$',default=x/(y^2 + x+1),width=60),
    P0=input_box(default=( n(.6,digits=3),n(.2,digits=3 )),label='$(x_0,y_0)$',width=35),
    parameter=input_box(default=(0,4),label='Interval for $a$',width=35),#,auto_update=False
    iterasyon=input_box(default=200,label='Quality Parameter#1',width=12),#number of iterations: 400 is good
    subint=input_box(default=200,label='Quality Parameter#2',width=12),#number of subinterval: 500 is good
    figs=slider(2, 20, 1, default=6,label='Figure size'),
    ar=slider([2^n for n in [-10..10]], default=1,label='Aspect Ratio' )
):
    paramin=parameter[0]
    paramax=parameter[1]
    kk=(paramax-paramin)/subint
    paralist=[paramin+i*kk for i in range(1+subint) ]
    biflistx=[]
    biflisty=[]
    for j in range(len(paralist)):
        nok=[(P0[0],P0[1] )]#list of points
        for i in range(iterasyon):
            nok.append((f(x=nok[i][0],y=nok[i][1],a=paralist[j] ),g(x=nok[i][0],y=nok[i][1],a=paralist[j])))
        nok_x=[ nok[-m][0] for m in [1..(iterasyon*2/5)] ]
        nok_y=[ nok[-m][1] for m in [1..(iterasyon*2/5)] ]
        for mm in range(len(nok_x)):
            biflistx.append( [paralist[j] ,nok_x[mm]])
            biflisty.append( [paralist[j] ,nok_y[mm]])
    biflistgrafx=list_plot(biflistx, plotjoined=False, color='grey',alpha=.3,pointsize=1)
    biflistgrafy=list_plot(biflisty, plotjoined=False, color='red',alpha=.3,pointsize=1)
    graf=biflistgrafx+biflistgrafy
    show(graf,frame=False,axes=True,axes_labels=['$a$','$x_n,y_n$'],figsize=figs,aspect_ratio=ar)

Interact: please open in CoCalc