mu = +1;
f = @(t,a) [ (mu - a(1).^2 - a(2).^2).*a(1) - a(2); (mu - a(1).^2 - a(2).^2).*a(2) + a(1)];
[t,a] = ode45(f,[0 100],[0.5 0.5]);     % Runge-Kutta 4th/5th order ODE solver
plot(a(:,1),a(:,2));
