% resevanje Van del Polove enacbe
% y'' = mu * (1-y^2) * y' - y                  (1)
% mu je parameter, ki doloca dusenje

% enacbo 2. reda prevedemo na sistem enacb 1. reda
% uvedemo z:=y' in Y := [y; z] in upostevamo (1)
%
% Y' = [y'; z']
%    = [y';   mu * (1-y^2)    * y'   - y]
%    = [Y(2); mu * (1-Y(1)^2) * Y(2) - Y(1)]

mu = 1;
vdp = @(t,Y) [Y(2); mu * (1-Y(1)^2) * Y(2)-Y(1)];
[t,Y] = ode45(vdp,[0,20],[2;0]);

% narisemo resitev
plot(t,Y(:,1),'-ro',t,Y(:,2),'-bo')
title(['resitev Van der Polove enacbe (mu=',num2str(mu),') z ode45']);
xlabel('t')
ylabel('resitev Y')
legend('y','z')