function ds = demo_odesoustava()
%
% Porovnani Runge-Kuttovy metody ctvrteho radu a Eulerovy metody:
% Reseni soustavy dvou rovnic druheho radu - problemu dvou teles
%

%%  Konecny cas a krok h
T = 20;  

h = 0.05;
%h = 0.01;
%h = 0.005;

%%  Pocatecni podminky:  u = [x  dx   y  dy ]
u0 = [1 0 0 1]; 

%%  Definice funkce f, plati ze u'=f(t,u)
f=inline('[w(2),-w(1)*(w(1)^2 + w(3)^2)^(-3/2),w(4), -w(3)*(w(1)^2 + w(3)^2)^(-3/2)]','w');

%%  Celkovy pocet kroku
pocetKroku = T / h; % celkovy pocet kroku

%%  Inicializace obrazku
figure;
axis([-1.5 1.5 -1.5 1.5]);
axis square
xlabel('X');
ylabel('Y');

%%  Startujeme s pocatecni podminkou
uR = u0; % RK4
uE = u0; % Euler

for i=0:pocetKroku
   
   disp(['krok ' num2str(i) ' z ' num2str(pocetKroku)])
   
   %%  Spocteme cas (t_n) odpovidajici aktualnimu kroku 
   t = i/pocetKroku * T; 
   pause(0.01);
   
   %%  Zobrazime aktualni polohu
   hold on;
   plot(uR(1),uR(3),'r.'); 
   plot(uE(1),uE(3),'b.'); 
   hold off;
   
   %%  Jeden krok Runge-Kutta 4
   k1 = f(uR);
   k2 = f(uR+h/2*k1);
   k3 = f(uR+h/2*k2);
   k4 = f(uR+h*k3);
   uR = uR + h/6 * ( k1 + 2*k2 + 2*k3 + k4);
   
   %%  Jeden krok Euler
   k1 = f(uE);
   uE = uE + h * k1;

   
end;

end
