glucose.mws

 >

Hopf bifurcation in the Glucosys model given by

,

 > restart; field := [-x+a*y+x^2*y, b-a*y-x^2*y];

 >

 > xs := b; ys := b/(a+b^2);

 > with(DEtools):with(plots): go1 := proc(pt) global ode, pp: ode := diff(x(t),t)=eval(field[1], [x=x(t), y=y(t)]), diff(y(t),t)=eval(field[2], [x=x(t), y=y(t)]): pp := dsolve({ode, x(0)=pt[1], y(0)=pt[2]},numeric); odeplot(pp, [x(t), y(t)], 150..180, view=[0..x1, 0..y1], thickness=5, numpoints=300, color=blue): end:

 >

 > go := proc() global a,b; local pic1, pic2; pic1 := DEplot([diff(x(t),t)=eval(field[1], [x=x(t), y=y(t)]), diff(y(t),t)=eval(field[2], [x=x(t), y=y(t)])], [x(t),y(t)], t=0..40,x=0..x1,y= 0..y1, map(e->[x(0)=e[1], y(0)=e[2]], pts), title=sprintf("a = %a;   b = %a",a,b) ,linecolor=1,arrows=MEDIUM, stepsize=0.04, thickness=2): pic2 := NULL: display(go1([1,1]),pic1 ,plot([b/(a+x^2),x/(a+x^2)], x=0..x1, color=black, linestyle=3, thickness=3) ): end:

 > a := 'a': b:='b':

 > pts := [[xs, ys+0.19], [1, 0], [0,2], [2,3.3]]; x1 := 3; y1 := 3;

 >

 > a := 0.1: b := 0.3:go();

 > a := 0.1: b := 0.4:go();

 > a := 0.1: b := 0.5:go();

 >

 > pic:=NULL: for b from 0.2 to 1 by 0.025 do go(); pic := pic, %: end:

 >

 > display(pic, insequence=true);

 >

 >