% Resavamo jednostavan Kosijev problem x''+x=0, x(0)=0, x'(0)=1, cije je jedinsteno
% resenje funkcija x(t)=sin(t). Resavamo tako sto prvo prebacimo u sistem:
% X' = [x1', x2'] = [x2, -x1], X(0) = [0, 1]. Kodiramo RK4.
clf

mode = 1; % 0 za crtanje grafika funkcije, 1 za crtanje fazne trajektorije

upper = 10;
t0 = 0; x0 = [0 1]';
h = 0.1; 
t = t0; 
x_osa = x0; % ovo je lista vektora vrednosti X(t) = [x1(t), x2(t)] = [x(t), x'(t)] koji trazimo
y1 = x0(1, 1);
y2 = x0(2, 1);
t_osa = t0:h:upper;
T = length(t_osa);

while (length(x_osa) < T) % RK4 metoda
    k1 = h * [y2 -y1]';
    k2 = h * ([y2 -y1]' + k1/2);
    k3 = h * ([y2 -y1]' + k2/2);
    k4 = h * ([y2 -y1]' + k3);
    
    k = [y1 y2]' + (k1 + 2*k2 + 2*k3 + k4)/6;
    y1 = k(1, 1);
    y2 = k(2, 1);
    x_osa = [x_osa [y1 y2]'];
    t = t + h;
end

x_osaA = sin(t_osa); % analiticko resenje

switch mode
    
    case 0

        plot(t_osa, x_osa(1, :)); hold on % plotujemo samo prvu vrstu, jer je x1(t)=x(t) nasa trazena funkcija
        plot(t_osa, x_osaA); hold on
        legend('RK4', 'Analiticko', 'Location', 'northwest');

    case 1
        
        plot(x_osa(1, :), x_osa(2, :)); hold on
        plot(sin(t_osa), cos(t_osa)); hold on % analiticko resenje za fazni portret. Dobija se kruznica parametrizovana sa (sin(t), cos(t))
        legend('RK4', 'Analiticko', 'Location', 'northwest');

end