%% Harj. 2, teht. 15 Doku
type diffsyst
clear  % Joskus on ihan "sallittua" Juhan varoituksista
       % huolimatta. :-)
       % (Varoitukset ehk? koskevat enemm?n Juhan massa-
       % ajoja teht?vien tarkistuksessa.)
%%
% Alkuarvot:
xyz0=[-10; -4; 12];

%% Ratkaistaan:
%  Muutellaan alkuarvoja lis??m?ll? 0.1 alkup. vektoriin
%  kullakin kierroksella. 
%  Kolmen sarakkeen tulosvektorit on k?tev?? 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??n 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
%% Tehd??n viel? animaatio.
figure
plot3(XYZ(:,1),XYZ(:,2),XYZ(:,3));
box;grid on