# Question:How to plot a prettier bifurcation diagram using Maple?

## Question:How to plot a prettier bifurcation diagram using Maple?

Maple

Dear All,

The following code plots the bifurcation diagram for a three-dimensional continuous dynamical system as a variable Re varies. However, the resulting plot (by pointplot command) is rather ugly, comparing with other bifurcation diagrams, see attached. Could anyone point me out how to improve its looking?

 > restart: with(plots): with(DEtools): with(plottools):with(LinearAlgebra): doSol:=proc( R )              local t_start:=2500,                    t_end:=5000,                    b:=-3,                    c:=3,                    v1:=1,                    f:=-4,                    v2:=2.0,                    omega:=0.1*sqrt(R),                    epsilon:=evalf(1/R),                    k:=0.1,                    kH:=(c+1)/(b-1),                    sys, evs, w, M, T, i, tt, solA, N_alpha,                    r_center, z_center, Zshift, alpha, alpha1,                    alpha2, del_alpha, m, Z, Rshift, RR;              sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u1(t),                   diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u2(t),                   diff(z(t),t)=z(t)*(kH*v1+c*u1(t)^2+c*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):              solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},                           {u1(t),u2(t),z(t)},                           numeric, method=rkf45, maxfun=10^7,                           events=[[[u1(t)=0, u2(t)>0], halt]]                          );              evs:=Array():              evs(1,1..4):=Array([t,u1(t),u2(t),z(t)]);              interface(warnlevel=0):              for i from 2 do                  w:=solA(t_end):                  if   rhs(w[1])=t_start                  then break;                  end if;              end do:              RR:=M[i..,3]: Z:=M[i..,4]:              N_alpha:=numelems(RR):              r_center:=add(RR[i],i=1..N_alpha)/N_alpha:              z_center:=add(Z[i],i=1..N_alpha)/N_alpha:              Zshift:=Z-~r_center: Rshift:=RR-~z_center:              alpha:=(arctan~(Zshift/~Rshift)+(Pi/2)*sign~(Rshift))/~(2*Pi)+~0.5:              return alpha;   end proc:
 > i_start:=125: i_end:=250: i_delta:=0.1: M:=[seq( [j, doSol(j)], j=i_start..i_end, i_delta)]:
 > S:=seq([Vector(numelems(M[i,2]),fill=M[i,1]),M[i,2]], i=1..numelems(M),1): display(seq( pointplot(S[i],              symbolsize=1,              symbol=point) ,              i=1..numelems(M)),              axes=boxed,              view=[i_start..i_end,0..1],              size=[650,400],              axesfont= ["TimesNewRoman", 16],              labels=["Re",phi[n]],              labelfont = ["TimesNewRoman", 16],              labeldirections=[horizontal, vertical]);

Bifurcation_diagram.mw

Thank you.

Very kind wishes,

Wang Zhe

﻿