function doe = demo_odesoustava_ellipt()
%
% Demonstrace Runge-Kuttovy metody ctvrteho radu pro reseni
% soustavy dvou rovnic druheho radu - problemu dvou teles
%

T = 10 ;  % budeme integrovat do 10 sekund

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

figure;
axis([-0.5 1.5  -1 1]);
axis square
hold on;
plot([-0.5 1.5],[0 0],'k--');
plot([0 0],[-1 1],'k--');
xlabel('X');
ylabel('Y');
hold off;

%%  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');

t = 0;
i = 0;
while(t<T)

    %%  vzdalenost od stredu
    r = (u(1)^2+u(3)^2)^(1/2);

    %%  delka casoveho kroku
    h = 1e-2;      %  pevny
    % h = 5e-3;      %  pevny 
    % h = 1e-3;      %  pevny
    % h = 1e-1*r^2;  %  "adaptivni"
   
   disp(['krok ' num2str(i) ': t=' num2str(t) ', r=' num2str(r) ', h=' num2str(h) ]);
   
   %% spocteme cas (t_n) odpovidajici kroku 
   i = i+1;
   t = t+h;
   pause(0.01);
   
   %% zobrazime aktualni polohu
   hold on;
   plot(u(1),u(3),'r.'); 
   hold off;
   
   %% jeden krok Runge-Kutta 4
   k1 = f(u);
   k2 = f(u+h/2*k1);
   k3 = f(u+h/2*k2);
   k4 = f(u+h*k3);
   u = u + h/6 * ( k1 + 2*k2 + 2*k3 + k4);
  
   
end;


end
