% ---------------------------------------------------------------------------
   % --------- ODE1 ------------------------------------------------------------ 
   % Standardni solver za resavanje ODJ primenom Euler-ove meode (Ojlerove metode). 
   % Poziv funkcije:
   % y_out = ODE1(F,t_0,h,t_final,y_0)  
   % Ojlerovom metodom se, sa fiksiranim korakom h na intervalu  t0 <= t <= t_final,
   % resava Kosijev zadatak:
   %      dy/dt = F(t,y)
   %      y(t0) = y0
   % Rezultat pozivanja funkcije je vektor y_out
   % ---------------------------------------------------------------------------

    function y_out = ode1(F,t0,h,tfinal,y0) 
      y = y0;
      y_out = y;
      for t = t_0 : h : t_final-h
         k = F(t,y);
         y = y + h*k;
         y_out = [y_out; y];
      end
    end

  % --------------------------------------------------------------------------------
  % Zadatak 1
  % Posmatrajmo Kosijev zadatak
  % dy/dt = 2*y; 
  % y(0) = 10
  % 0 <= t <= 3
  % --------------------------------------------------------------------------------
  % >> F = @(t,y) 2*y;
  % >> y0 = 10;
  % >> t_0 = 0;
  % >> t_final = 3;
  % >> h = 1;      
  % >> y = ode1(F,t_0,t_final, y0)
  % ako zelimo da animiramo sva tri koraka metode, kucamo naredbu
  % >> euler_steps                   
  % --------------------------------------------------------------------------------  
  % Zadatak 2 
  % Slozeni kamatni racun
  % >> p = 0.06; 
  % >> F = @(t,y) p*y;
  % >> t_0 = 0;
  % >> h = 1/12;
  % >> t_final = 10; 		   % 10 godina
  % >> y0 = 1000;                    % 1000 rsd (eur, dolar, ...)
  % >> y = ode1(F,t_0,t_final, y0)
  
  % format bank    % novcane jedinice su dinari i pare, evri i centi, ... (2 decimale). 
  % grafik + poredjenje sa prostim kamatnim racunom y1 = y0 + y0+p*t;
  % >> t = (0:h:t_final)';  
  % >> plot(t,[y y1],'.')         
  % --------------------------------------------------------------------------------


   % --------------------------------------------------------------------------------
   % ---------------------- ODE2 ---------------------------------------------------- 
   % Standardni solver za resavanje ODJ primenom meode pravougaonika
   % Poziv funkcije:
   % y_out = ODE2(F,t_0,h,t_final,y_0)  
   %  metodom sa fiksiranim korakom h na intervalu  t0 <= t <= t_final
   % resava sledeci Kosijev zadatak:
   %      dy/dt = F(t,y)
   %      y(t0) = y0
   % ---------------------------------------------------------------------------

    function y_out = ode2(F,t0,h,tfinal,y0) 
      y = y0;
      y_out = y;
      for t = t_0 : h : t_final-h
         k1 = F(t,y);
         k2 = F(t+h/2, y+h*k1/2);
         y = y + h*k2;
         y_out = [y_out; y];
      end
     end


  % Zadatak 3:
  % Posmatrajmo Kosijev zadatak kojim se opisuje sirenje plamena sibice
  % dy/dt = sqrt(1-y^2);     
  % y(0) = y0
  % 0 <= t <= 2/t_0
  % Tacno resenje Kosijevog problema je y(t) = sint;
  % --------------------------------------------------------------------------------
  % podaci.m
  % >> F = @(t,y) sqrt(1-y.^2);
  % >> y_0 = 0;
  % >> t_0 = 0;
  % >> t_final = pi/2;
  % >> h = pi/32;     
  % y = ode2(F,t_0,t_final, y_0)
  % t = (0:h:t_final)';
  % plot(t,y,'o-')
  % set(gca, 'xtick', [0:pi/8:pi/2], 'xticklabel', {'0','\pi/8', '\pi/4', '3\pi/8', 'pi/2'}) 
  % axis([0 pi/2 0 1]) 
  % xlabel('t'), ylabel('y')       


  % ako pogledamo vrednosti vektora y videcemo koliko resenja ove metode odstupaju od tacnog resenja 
  % buduci da znamo vrednosti funkcije sin(t).



  % Zadatak 4 (isti kao 2):
  % >> F = @(t,y) 2*;
  % >> t0 = 0;
  % >> h = 1;
  % >> tfinal = 3;
  % >> y0 = 10;
  % >> y = ode2(F,t_0,t_final, y_0)
  % >> midpoint_steps  
 
  % VEZBA : POGLEDATI METOD ode2t !!!

 

   % --------- ODE4 ------------------------------------------------------------ 
   % Standardni solver za resavanje ODJ primenom meode Runge-Kutta.
   % Poziv funkcije:
   % y_out = ODE4(F,t_0,h,t_final,y_0)  
   % Runge-Kutta metodom sa fiksiranim korakom h na intervalu  t0 <= t <= t_final
   % resava sledeci Kosijev zadatak:
   %      dy/dt = F(t,y)
   %      y(t0) = y0
   % ---------------------------------------------------------------------------

    function y_out = ode4(F,t0,h,tfinal,y0) 
      y = y0;
      y_out = y;
      for t = t_0 : h : t_final-h
         k1 = F(t,y);
         k2 = F(t+h/2, y+h*k1/2);
         k3 = F(t+h/2, y+h*k2/2);
         k4 = F(t+h, y+h*k3);
         y = y + h*(k1 + 2*k2 + 2*k3 + k4)/6;
         y_out = [y_out; y];
      end
     end


  % Zadatak 1:
  % Posmatrajmo Kosijev zadatak kojim se opisuje sirenje plamena sibice
  % dy/dt = y^2 - y^3
  % y(0) = y_0
  % 0 <= t <= 2/t_0

  % --------------------------------------------------------------------------------
  % podaci.m
  % F = @(t,y) y.^2-y.^3;
  % y_0 = 0.01;
  % t_0 = 0;
  % t_final = 2/y_0;
  % h = t_final/500;      % 500 koraka
  % y = ode4(F,t_0,t_final, y_0)
  % t = (0:h:t_final)';
  % plot(t,y,'.')         
  % primeticemo da plan pocinje da raste postepeno, zatim se naglo rasiri, 
  % nakon cega ostaje fiksiranog poluprecnika

  % resiti Zadatak 1 sa drugacijim ulaznim podacima:
  % t_0 = 0;
  % y_0 = 1/1000;
  % t_final = 2/y_0;
  % h = t_final/2000;
  % za koju vrednost parametra t plamen dostize poluprecnik 0.9?
  % ---------------------------------------------------------------------------------
 

  % Zadatak 2:
  % --------------------------------------------------------------------------------
  % podaci.m
  % F = @(t,y) 1 + y.^2;
  % y_0 = 0;
  % t_0 = 0;
  % t_final = 2;
  % h = t_final/500;      % 500 koraka


  % STEPEN TACNOSTI METODA !
  % stepen tacnosti mozemo odrediti tako sto izracunamo vrednost problema sa korakom h, zatim sa korakom h/2 i uporedimo sa tacnim resenjem

  % Posmatramo sledeci problem:
  % Integral_0^1 (1/(1+t)^2)dt
  % Tacno resenje integrala je 1/2.

  function p = orderx(odex)
    
    vexact = 0.5;
    F = @(t,y) 1/(1+t)^2;
    t0 = 0;
    tfinal = 1;
    y0 = 0;
    h = 0.1
    y_out = odex(F,t0,h,tfinal,y0);
    v1 = y_out(end)
    h = h/2
    y_out = odex(F,t0,h,tfinal,y0);
    v2 = y_out(end)
    ratio = (v1 - vexact)/(v2 - vexact)
    p = round(log2(ratio));
   
 end
   
  % primer: 
  % orderx(@ode1)
  % orderx(@ode2)
  % orderx(@ode4)



 % ---------- ode23 ----------------------
 % umesto da definisemo korak, mozemo da definisemo tacnost koju zelimo da postignemo.
 % Bogacki - Shampine odred 3(2) 
 % Koriste 3 vrednosti funkcije kako bi postigli tacnost 2.
 % k1 = f(t,y)
 % k2 = f(t+h/2, y+h/2*k1);
 % k3 = f(t+3/4h, y+3/4*h*k2);
 
 % y = y + h/9*(2*k1 + 3*k2 + 4*k3);
 
 % k4 = f(t+h,y);
 % e = h/72*(-5*k1+6*k2+8*k3-9*k4); 

 % >> F = @(t,y) y;
 % >> ode23(F,[0 1],1)     % zanima nas resenje problema dy/dt=y na segmentu [0,1] ako je y0=1;
                           % nismo napisali ni jedan argument za izlaz
 		           % samo ce nacrtati resenje
 % >> [t,y] = ode23(F,[0,1],1)
		           % vratice resenje ako stavimo izlazne argumente
                           % sam bira korak za t.

 % tolerancija za resenje je 10^(-3).

 % Primer:
 % a = 0.25;
 % y0 = 15.9
 % F = @(t,y) 2*(a-t)*y^2;
 % [t,y] = ode23(F, [0,1],y0);
 % plot(t,y,'-',t,y,'.','markersize', 0.24)
 % h = diff(t);



 % ----------- ode45 -----------------------
 % Dormand - Prince su reda 5(4)
 % >> F = @(t,y) y;
 % >> tspan = (0:0.1: 1)'
 % >> [t y] = ode45(F,tspan,1)
 % >> plot(t,y,'o-') 
 % ako zelimo da proverimo tacnost
 % format short e
 % >> exp(t)-y         % tacno resenje
 % ans = ......
 % a = 0.25;
 % y0 = 15.9;
 % F =@(t,y) 2*(a-y)*y^2;
 % [t y] = ode45(F,[0,1],y0)     % ode45 sam bira svoj korak
 %