% metody typu Runge-Kutta

x0 = 0;
xE = 2;
Y0 = [1;1];
h  = 0.001;

% presne reseni
[x1,y1] = RK_method(x0,xE,Y0, 0.0001,1);
% RK1
[x2,y2] = RK_method(x0,xE,Y0, 1/250,1);
% RK2
[x3,y3] = RK_method(x0,xE,Y0, 1/250,2);
% RK3
[x4,y4] = RK_method(x0,xE,Y0, 1/250,3);


% Zobraz vysledky
figure (1);
hold on;
plot(x1,y1(1,:),'-k');
plot(x2,y2(1,:),'or');
plot(x3,y3(1,:),'ob');
plot(x4,y4(1,:),'og');

xlabel('x_n',"fontsize", 18);
ylabel('y_n',"fontsize", 18);
h=get(gcf, "currentaxes");
set(h, "fontsize", 18, "linewidth", 1);
grid on;
hold off;

function F = F(x,Y)
    m = 0.1;
    k = 4.2;
    d = 0.01;
    F = [ Y(2);
          1/m * ( 1/(10*x+1)*sin(50*pi*x) - d*Y(2) - k*Y(1) )
         ];
end

function [x,Y] = RK_method(x0,xE,Y0,h,order)
    x = x0:h:xE;    % x-ove body aproximace
    N = length(x);  % pocet iteraci
    Y(:,1) = Y0;    % pocatecni podminka

    if      (order == 1) % Explicitni Euler (RK1)
        for i=1:N-1
            Y(:,i+1) = Y(:,i) + h * F(x(i),Y(:,i)); 
        end
    elseif  (order == 2) % Collatzova metoda (RK2)
        for i=1:N-1
            k1 = Y(:,i) + 0.5*h * F(x(i),Y(:,i));
            k2 = F(x(i) + 0.5*h, k1);            
            Y(:,i+1) = Y(:,i) + h * k2;
        end        
    elseif  (order == 3)
        for i=1:N-1
            k1 = F(x(i),Y(:,i));
            k2 = F(x(i)+0.5*h, Y(:,i) + 0.5*h*k1);
            k3 = F(x(i)+h, Y(:,i) + 2*h*k2 -h*k1);
            Y(:,i+1) = Y(:,i) + h/6*(k1 + 4*k2 + k3);
        end
    end

end