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

omega = 1;
a     = 1;
b     = 2;
k     = 10;

x0 = 1;
y0 = 1;


p = [omega  a  b  k];
v0 = [x0 y0];

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

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

figure(1)
plot(t,v);

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(@flinsys,[0 T],v0,opts,p);

figure(2)
plot(t,v)

data = [t v omega*t];

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