function f = vanderpol_oscillator(x0,mu)

T = 3*pi;
%dt = 0.05;
%t = 0:dt:T;

OPT = odeset('RelTol',1e-3,'AbsTol',[1e-3]);

[t,x] = ode45(@vdp_oscillator,[0 T],x0,OPT,mu);


plot(0,0,'*')
set(gca,'XLim',[-5 5],'Ylim',[-5 5], ...
    'DataAspectRatio',[1 1 1])
hold on
for i = 1:length(t)
    plot(x(1:i,1),x(1:i,2),'b-','LineWidth',2);
    pause(0.01)
end
xlabel('x(1)')
ylabel('x(2)')


set(gca,'FontSize',16,'FontWeight','b','FontAngle','i')
if mu == 0
    title('Harmonic Oscillator');
else
    title('van der Pol Oscillator');
end
xlabel('x(1)')
ylabel('x(2)')

f = 'EXIT_SUCCESS';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dx] = vdp_oscillator(t,x,mu)

dx = zeros(2,1);

dx(1) = x(2);
dx(2) = -x(1) - mu * (x(1)^2-1)*x(2);