function f = linearized_vanderpol(x0,mu)

T = 3*pi;
dt = 0.1;
t = 0:dt:T;

Kolor = rand(1,3);

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

[t,x] = ode45(@vdp_oscillator,t,x0,OPT,mu);
[t,z] = ode45(@linear_vdp_oscillator,t,x0,OPT,mu);

subplot(121)

plot(0,0,'*')
set(gca,'XLim',[-4 4],'Ylim',[-4 4], ...
    'DataAspectRatio',[1 1 1])
hold on

for i = 1:length(t)
    plot(x(1:i,1),x(1:i,2),'Color',Kolor,'LineWidth',2);
    pause(0.01)
end
set(gca,'FontSize',12,'FontWeight','b','FontAngle','i')

if mu == 0
    title('Harmonic Oscillator');
else
    title('van der Pol Oscillator');
end

xlabel('x(1)')
ylabel('x(2)')





subplot(122)
plot(0,0,'*')
set(gca,'XLim',[-4 4],'Ylim',[-4 4], ...    
    'DataAspectRatio',[1 1 1])
hold on
for i = 1:length(t)
    plot(z(1:i,1),z(1:i,2),'Color',Kolor,'LineWidth',2);
    pause(0.01)
end

set(gca,'FontSize',12,'FontWeight','b','FontAngle','i')
if mu == 0
    title('Harmonic Oscillator');
else
    title('Linearized van der Pol Oscillator');
end
xlabel('z(1)')
ylabel('z(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);

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

dx = zeros(2,1);

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