%%
%% Listing 7.10 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

% oscylator Van der Pola

% parametry zadania
epsilon = 0.8;
Tmin = 0; 
Tmax = 60; 
Y0 = [0.01; 0]; % warunek początkowy

% t - zestaw $N$ punktów, w których chcemy wyznaczyć rozwiązanie y na przedziale [Tmin, Tmax]
N = 1000;
t = linspace (Tmin, Tmax, N);

% rozwiąż! Y - tablica wartości rozwiązania w punktach t

if exist('OCTAVE_VERSION') 
% wersja dla Octave
	vanderpol = @(Y,t) vanderpoleps(t,Y,epsilon);
	[Y, info, msg] = lsode(vanderpol, Y0, t);
	if (info ~= 2)
		disp(['Problemy w dzialaniu lsode: ', msg]);
	end
else
% wersja dla MATLAB-a
 	vanderpol = @(t,Y) vanderpoleps(t,Y,epsilon);
	[t, Y] = ode45(vanderpol, t, Y0);
end

% wizualizacja i zapis wykresów do pliku
filename = ['vanderpol', 'T-', num2str(Tmax,'%g'), ...
		'epsilon-', num2str(epsilon,'%g')]; % główny człon nazwy pliku

plot(t, Y(:,1)); % wykres rozwiązania $y(t)$
title('Oscylator Van der Pola');
xlabel('czas (t)'); 
legend('y(t)');
print('-depsc', [filename, '-sol.eps']);

input('Nacisnij dowolny klawisz...'); % czekaj na naciśnięcie klawisza!

plot(Y(:,1), Y(:,2)); % wykres fazowy rozwiązania $(y(t),y'(t))$
title('Krzywa fazowa oscylatora Van der Pola');
xlabel('y_1(t)'); ylabel('y_2(t)');
print('-depsc', [filename, '-phase.eps']);
