% Simple plot of forced, damped Duffing equation
%
% d^2 y/dt^2 = y-y^3-delta*dy/dt + gamma*cos(om*t);
%
% Left click for initial conditions, Runs until beep sounds. 
% Right click to stop or left click for another orbit.
% Need function  fduffeq.m

clear
clf
hold on
h = animatedline('LineStyle','none','Marker','.','MarkerSize',2);
axis([-2 2 -2 2]);
om=1.0;
per=2*pi/om;
t1=0;
[x0 y0 bb]=ginput(1);  
addpoints(h,x0,y0);
drawnow update;
x0=-1;
% for j=1:1
%   y0=-1+(j-1)*2/20;
   while bb<2, 
      for k=1:1000
       t2=t1+per;
       [T,Y] = ode45('fduffeq',[t1 t2],[x0 y0]); 
       [nrow,ncol]=size(T);
       x0=Y(nrow:nrow,1);
       y0=Y(nrow:nrow,2);
       % plot(x0,y0), 
       addpoints(h,x0,y0);
       drawnow update;
      
       t1=t2;
      end
    t1=0;
    beep
    [x0 y0,bb]=ginput(1);
  end
hold off;
