opts = odeset('reltol',1e-6,'abstol',1e-8);

omega = 1;
a     = 5;
e     = 1/10;

x0 = 1;
y0 = 1;


p = [omega  a  e];
v0 = [x0 y0];

T = 2*pi/omega;
Ttrans = 80*T;

% Transient...
[t,v] = ode45(@wnsys,[0 Ttrans],v0,opts,p);

figure(1)
plot(t,v);
grid on

v0 = v(end,:);

% Use the last data point of the previous solution
% as the initial condition, and solve for one period.
% The should be a good approximation to the periodic orbit.
[t,v] = ode45(@wnsys,[0 T],v0,opts,p);

figure(2)
plot(t,v)
grid on

data = [t v omega*t];

save('wn.dat','data','-ascii')