%% Harj. 2, teht. 15 Doku
% Kirjoitetaan differentiaaliyhtälösysteemin koodi tiedostoon diffsyst.m
%{
function dy  = diffsyst( t,x )
dy=[    -10*x(1)+28*x(2);       ...
        10*x(1)-x(2)-x(1)*x(3); ...
        -8/3*x(3)+x(1)*x(2)     ...
   ];
%}
% Huomaa, että muuttuja t tarvitaan, vaikka sitä ei tässä (autonomisen 
% yhtälön tapauksessa käytetä mihinkään. (Kutsusyntaksin oltava sama, jotta ...
    kutsuva ode-rutiini voi toimia kaikenlaisille yhtälöille.)
%%
% Alkuarvot:
xyz0=[-10; -4; 12];

%% Ratkaistaan:
%  Muutellaan alkuarvoja lis"a"am"all"a 0.1 alkup. vektoriin
%  kullakin kierroksella. 
%  Kolmen sarakkeen tulosvektorit on k"atev"a"a tallettaa
%  solutaulukkoon ("cell array"). Se voidaan ajatella 
%  indeksoiduksi muuttujaksi, jota indeksoidaan 
%  aaltosuluilla.
N=2;

for k=0:N
  xyz0=(1+k*0.1)*ones(size(xyz0));
 [T,XYZ]=ode45(@diffsyst,0:.01:40,xyz0);
 XYZM{k+1}=XYZ;
end 
%% Verrataan pienen muutoksen vaikutusta ratkaisuun.
max(abs(XYZM{1}(:,1)))
max(abs(XYZM{1}(:,1)-XYZM{2}(:,1)))
subplot(2,1,1)
  plot(T,XYZM{1}(:,1))
subplot(2,1,2)
  plot(T,XYZM{2}(:,1))
shg  
%%  
figure
plot(T,XYZM{2}(:,1)-XYZM{1}(:,1))
shg
%% Piirret"a"an kaikki kuvat viimeksi saaduilla arvoilla.
figure
subplot(3,1,1)
plot(T,XYZ(:,1))
subplot(3,1,2)
plot(T,XYZ(:,2))
subplot(3,1,3)
plot(T,XYZ(:,3))
shg
%% 
figure
plot3(XYZ(:,1),XYZ(:,2),XYZ(:,3));
box;grid on
%% Animaatio komennolla lorenz
% Kokeile, lue doc lorenz ja Molerin NCM-kirjan sivuilta 
% 202-203 "Lorenz attractor"
% 