% define parameters
% h= timestep , a = driving amplitude , b = natural parameter
% c = damping parameter , w = driving frequency
%a=0.9     % periodic
%a=1.07    % period doubling
%a=1.15    % chaos
%a=1.35    % periodic
%a=1.45    % period doubling
%a=1.47    % period doubling
a=1.50;     % chaos
b=1.0;
c=0.5;
w=0.66666666;
h=0.01;
% define length of main loop
num=100000;
% arbitrary starting point
x=0.0;
y=0.0;
% open a window
figure('Position',[50,50,600,600]);
% set up plot
p=plot(x,y,'EraseMode','none', 'Color', 'red');
axis([-pi,pi,-4.0,4.0]);
xlabel('Position');
ylabel('Velocity');
title('Non-Linear Oscillator - Phasespace Plot');
% use loop to generate successive points in time
% and plot points in real time
for j=1:num
  % Runge-Kutta method
  f1=y;
  g1=osc1(j*h,x,y,a,b,c,w);
  f2=y+h*g1/2.0;
  g2=osc1((j+0.5)*h,x+h*f1/2.0,y+h*g1/2.0,a,b,c,w);
  f3=y+h*g2/2.0;
  g3=osc1((j+0.5)*h,x+h*f2/2.0,y+h*g2/2.0,a,b,c,w);
  f4=y+h*g3;
  g4=osc1((j+1)*h,x+h*f3,y+h*g3,a,b,c,w);
  x=x+h*(f1+2.0*f2+2.0*f3+f4)/6.0;
  % keep x = angular variable between -2*pi and +2*pi
  if (abs(x) > pi )
      x = x-2.0*pi*sign(x);
  end
  y=y+h*(g1+2.0*g2+2.0*g3+g4)/6.0;
  %get rid of transients before plotting
  if (j > 25000)
    set(p,'XData',x,'YData',y)
    drawnow
  end
end
