%%
%% Listing 7.11 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/
%%

% rozwiązanie równania Lorenza

sigma = 10;
r = 28;
b = 8/3;

% t - zestaw punktów, w których chcemy wyznaczyć rozwiązanie y 
% na przedziale [0,T]

T = 200;
N = 10000; 	% wybieramy bardzo dużo punktów, by mieć lepszą wizualizację skomplikowanego rozwiązania
t = linspace (0, T, N);

% y0 - warunek początkowy

y0 = [0.01; 0; 0];

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

if exist('OCTAVE_VERSION') 
% wersja dla Octave
	[y, info, msg] = lsode( @(y,t)(lorenz(t,y,sigma,r,b)), y0, t);
	msg
else
% wersja dla MATLABa
	[t, y] = ode15s( @(t,y)(lorenz(t,y,sigma,r,b)), y0, t);
end
% wizualizacja i zapis do pliku

plot(y(:,1), y(:,2), 'r-');
title('Rozwiazanie modelu Lorenza: rzut 2D');
xlabel('x(t)'); ylabel('y(t)');

print('-depsc', 'lorenz2d.eps');

plot3(y(:,1), y(:,2), y(:,3), 'b-');
title('Rozwiazanie modelu Lorenza: rzut 3D');
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)');

print('-depsc', 'lorenz3d.eps');
