@dharr it is my idea i don't know it is apply or not, but i think with your coding it will be easy, if we change all the equalibriom points to (1,0) and (-1,0) and in some of them both are exist (-1,#) it will be more easier for jacobian and phase portrait, last time you write a good programing for calculate that ,i used for finding the parameter i think it will be usfull i want used for all equalibriom point in once, and if this apply for each equalibriom  then when we call equalibriom points in phase portrait it will be clear which each equalibriom point what behavior have it, i did a paper before but they are not like this really i need help about that to find the equalibriom and phase portrait, and the your jacobian is great i liked really , this is one one of my old paper i did like that  the programing already in top but for finding parameter i mention another maple file, if the equalibriom are the simple is not efect the behavior so we need use the simple one  and 1,-1 are simplest one
and sorry for top i see that the picture are small i fixed that too

Findin-S-parameter=1.mw

@dharr i think you are write about the condition when he mention i think is wrong becuase the condition of equality 

diff(P,Omega) = -diff(Q,y), is not satisfy while i am trying to  search for more equation i found just one another paper which same other solved like this shape is something new  and this type don't have much paper which i am intresting to do it better than him, i try to solve another paper but i  didn't got same ode like he did but i did exact imaginary part condition when beta=-2*p at the paper, by my openion this  beta=0 as result say but i did same as paper did so in finding the conserved quantity  susccesfull and the condition satisfy but for phase portrait and jacobian is a litle bit hard this time contain two function this is my first time to do a system contain two function by programing and thus point, also i think if i replace the variable  like 

y(eta)*(2*chi*Omega(eta)^2 + p)=s[1]*V*U^2+V*s[2]

when s[1]=2*chi and s[2]=p something like that it will be better even for finding equalibriom points, in below i did my best try for finding ode and the rest but  iam stuck ? even i try to work more for the first equation at there how we can find thus equation to make them be zero ? can we do anything like that?

bi-2.mw

F-ode-Fpaper.mw

bifur.mw

if need any thing please let me know, 

second paper

first paper

@acer @Carl Love is not working there is any other way for doing this plot? i want to plot a good look of my chaotic plot but i just have this code which most time not working?

@sand15  the persion do that he is expert , if you have any idea please do

@dharr  it is true but for phase portrait specially for chaotic this make a problem becuase it chose a complex number  and when i used complex number chaotic not appear anything 

i tested on that S[1] and S[2] in complex equation which totally satisfy thanks for code

true.mw

@mmcdara right now i am reach on this plotting, and i am notic this which for work better we have to do something is ok if we don't make 3D explor option we can use 2D explor option for outcome of 3D shape becuase they related, so i need thus two part in after explor option i need the system and parameter outcome,  for plotting, also in case if we want to convert for other language of programing How we can do that i want plot that shape of 2D when i construct by explor option?

i can't jump in your code they sensitive a little bit

pp.mw

@dharr  let say S[1]=x+y+z
i want find S[1]=1, what is x+y+z which x,y,z are not be zero , let say x=1,y=1, so z=-1 becuase parameter sometime are a lot i want something automatically and randomly find them if x,y,z are be rond it will be better 

i think this time is clear, in the paper he did like that is totally clear 

@dharr  i just change letter and i solve another problem is not like this but idea is same  you say mu[1] mu[2] any letter just find element iside them which make them 1 or -1 it is easier for showing and for phase portrait it is best idea  which be 1,

you can use the picture system in case the exact solution and result are clear.

about plot i don't know how plot imaginary it mean  must have another structure or what i don't have idea  i can't go to global bifurcation while i am stuck in local , 

@mmcdara i hope you be able to to explore option for phase portrait too , i have work on it and with your comunity maple i think i can do the best job ever possible  thank you 

@mmcdara  Wow, you just explode  it is so emazing  really i have to work on it to improve this work i have 1000 question in my mind, i saw the explore option for 3D too, i have the code for easy function in case needed mention, mr accr did for me a long time ago, but for this it is in my mind i say i will ask in future but you did it now and it is far amazing this is made my day, 

i agree with you i saw the rung kutta method they used for sensitivity but  i am not sure how and i am not involve that part yet but for sensity Mr.rubien did another way which yet i didn't have much time for work on that too i have to write a paper and part by part i am stuck  becuase this is my first work about dynamic system, 
additionally, i have a lot question  and i will ask after , please if you can do for 3D it will be something new  my future work and i don't think now one done this 

@dharr  about point 1: if you watch S[1] in fact contain p[1],p[2],p[3] and k[3] inside S[1] i want set all of p and k and other parameter inside S[1] also for S[2] which the result be 1, if you watch the picture i mentioned the different picture, if you watch case1: you will see in system we have mu[1] and mu[2] but inside of them we have other parameter which author set them up arrange them some how somekind of solve all parameter use  to find mu[1] and mu[1] in equalibriuom point i have to say is (1,0) instead of (sqrt(S[1]/S[2],0) aso i have to say what is parameter inside S[1] and S[2] for getting this (1,0) 

point2: if you watch the equalibriom point we have three one which two of them are complex i don't know why not shown this is my first time i write a paper so i am looking for information stil about that i know first integral plot are great  but just show the point (0,0) how about (i,0) and (-i,0) i try ask Ai he say we need global phase portrait which i don't have idea how  can show imaginary part of phase portrait  i didn't see that, in paper when they have this case they just say we have (0,00 but i want to know why other not include

please if don't undrestand mention i will explain more and if need example 

@sand15 i can't upload why i don't know but i take a screenshot

@mmcdara i don't beliave that which you done this, this is amazing but i have work on it and fully undrestand the coding also i have a lot new idea about getting new shape of thus chaotic, and the code is well runing untill here 

YOUR MATLAB EXAMPLE

is not give me outcome i waited for 5 minute at least but no outcoming 

i have to work on it and used for other example also new idea if you have a time i can give you a example with new idea 

i just got a new article from my friend which really outcome are amazing and i want to ask a new question about it but in here i think is better to present becuase your outcome are great and quite close let change the system and new using function 

@dharr  i don't know why the picture in here are so different from matlab, i try to use matlab but at there i didn't give xi and otehr parameter i just write the system it self with k[1] and k[2] with thus option

K1 = 1.0;      % Cubic nonlinearity coefficient
K2 = -1.0;     % Linear restoring force coefficient  
epsilon = 1.0;  % Forcing amplitude
theta = 1.2;   % Forcing frequency

% Time span
tspan = [0 200];  % Long integration time for chaotic behavior
dt = 0.01;
t = 0:dt:tspan(2);

% Initial conditions
u0 = 0.1;  % Initial position
v0 = 0.1;  % Initial velocity
initial_conditions = [u0; v0];
%% Define the system of ODEs
duffing_system = @(t, y) [
    y(2);                                    % du/dt = v
    -K1*y(1)^3 + K2*y(1) + epsilon*sin(theta*t)  % dv/dt
];

and i get the good result  

@mmcdara  thank you so much i will wait , i don't undrestand the first coment where i should add and i paste your code there but nothing happen, and when i change the name of color in cyan to cyan or gold to Gold is give me error so colors are still have same name