% Resavamo jednostavan Kosijev problem x'=t+x, x(0)=e-1, cije je jedinsteno
% resenje funkcija x(t)=e^(t+1)-t-1. Recimo da trazimo resenje na intervalu [0, upper]
clf

upper = 5;
t0 = 0; x0 = exp(1) - 1;
h = 0.001; % biramo korak
t = t0; % krecemo od pocetnog vremenskog trenutka t0 i ici cemo sve do upper
x_osa1 = x0; 
x_osa2 = x0; % uzmemo listu vrednosti x koje treba da dobijemo, sa jednim elementom x0
y1 = x0; 
y2 = x0; % y nam predstavlja novu numericku vrednost aproksimacije. Na pocetku je to tacna vrednost x0
t_osa = t0:h:upper;
T = length(t_osa);

while (length(x_osa1) < T) % Ojlerova metoda, koristimo y1 i x_osa1
    y1 = y1 + h * (t + y1);
    x_osa1 = [x_osa1 y1];
    t = t + h;
end

t = t0;

while (length(x_osa2) < T) % RK4 metoda, koristimo y2 i x_osa2
    k1 = h * (t + y2);
    k2 = h * (t + h/2 + y2 + k1/2);
    k3 = h * (t + h/2 + y2 + k2/2);
    k4 = h * (t + h + y2 + k3);
    
    y2 = y2 + (k1 + 2*k2 + 2*k3 + k4)/6;
    x_osa2 = [x_osa2 y2];
    t = t + h;
end

x_osaA = exp(t_osa+1)-t_osa-1; % analiticko resenje
plot(t_osa, x_osa1); hold on
plot(t_osa, x_osa2); hold on
plot(t_osa, x_osaA); hold on

legend('Ojler unapred', 'RK4', 'Analiticko', 'Location', 'northwest');