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

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_end
                 then evs(i,1..4):=Array(map(rhs, w));
                      solA(eventclear);
                 else break;
                 fi
             od:

             interface(warnlevel=3):
             M:=DeleteRow(convert(evs,matrix),1):

             T:=M[..,1]:

             for i from 1 do
                 tt:=T[i]:
                 if   tt>=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]);