function dso = demo_simpleode()
%
% Ukazka nejjednodussich Runge-Kuttovych metod: Eulerovy, midpoint, Heunovy
%
%
% Resena dif. rovnice je  u'(t)=f(t,u)
% Zadavame  
%   f, tedy pravou stranu
%   presne reseni - pouze pro vykresleni
%   pocatecni podminku
%   konecny cas
%   delku kroku h

uloha = 1;

if (uloha==1)
    % Bakterie.
    f      = inline('(1+cos(t)).*y;','t','y'); 
    fexact = inline('exp(t + sin(t));','t'); 
    y0     = 1;
    T      = 5;
    h  = 0.1; 
    %h  = 0.05; 
    %h  = 0.01; 

elseif (uloha==2)
    % Obtiznejsi uloha.
    % Hranice pouzitelneho h pro Midpoint a Heuna je mezi 0.13 a 0.14.
    f      = inline('-15.*y','t', 'y'); 
    fexact = inline('exp(-15.*t)','t'); 
    y0     = 1;
    T      = 5;
    h  = 0.14; 
    %h  = 0.13; 
    %h  = 0.12; 
    
end

%%  Pole obsahujici t_i a presne reseni y_i
tA    = 0:h:T;
exact = fexact(tA);
nbodu = size(tA,2);

%%  Inicializujeme reseni jednotlivymi metodami (aplikace poc. podm.)
uE = y0;
uM = y0;
uH = y0;

%%  Zobrazime presne reseni a nastavime grafy aby se nemazaly,
%%  ale aby se jen pridavaly objekty do existujicich grafu
figure;
pex = plot(tA,exact,'k--'); % zobrazime presne reseni
hold on;

%%  Zobrazime pocatecni podminky
t = tA(1);
pE = plot(t,uE,'r+-');
pM = plot(t,uM,'g*-');
pH = plot(t,uH,'bd-');
legend({'presne','Euler','Midpoint','Heun'},'location','NW');

for i=1:nbodu-1,
    %%  V cyklu vzdy postoupime danou metodou o krok h
    %%  a nakreslime aktualni hodnotu y_i+1
    %%     (presneji receno spojene hodnoty y_i a y_i+1)
    
    t = tA(i);

    %%  Euler
    ui = uE;  % pomocna promenna pro vykresleni (spojeni s y_n)
    uE = uE + h*f(t,uE); 
    plot([t t+h],[ui uE],'r+-'); 

    %%  Midpoint
    ui = uM;
    uM = uM + h*f(t+h/2,uM+h/2*f(t,uM));
    plot([t t+h],[ui uM],'g*-'); 

    %%  Heun
    ui = uH;
    uH = uH + h/2*(f(t,uH)+f(t+h,uH+h*f(t,uH)));
    plot([t t+h],[ui uH],'bd-'); 

end;

hold off;

end

