Derivare numerica

function y = Euler(x0, a, n, f, y0)
% Functia care calculeaza solutia ecuatiei: y' = f(x, y) prin metoda Euler.
%   Date de intrare:
%       - x0 -> capatul stang al intervalului;
%       - a  -> lungimea intervalului;
%       - n  -> numarul de puncte;
%       - y0 -> conditia initiala;
%       - f  -> functia de integrat.
%   Date de iesire:
%       - y  -> vectorul aproximatiilor solutilor.
    y = zeros(n+1, 1);
    y(1) = y0;
    h = a / n;
    for i = 2 : n + 1
        x = x0 + (i - 1) * h;
        y( i) =y(i-1) + h * f(x, y( i-1 ) );
    endfor
endfunction


function y = PredictorCorector(a, b, n, y0, f)
%   Functia care calculeaza solutia ecuatiei
%y' = f(x, y) prin metoda Runge - Kutta
%   Date de intrare:
%       - a, b -> intervalul de integrare;
%       - n    -> numarul de puncte;
%       - y0   -> conditia initiala;
%       - f    -> functia de integrat y'=f(x,y).
% Date de iesire:
%       - y    -> vectorul aproximatiilor solutiei.
    y = Runge_Kutta(a, b, n, y0, f);
    h = (b-a) / n;
    x0=a;
    x1 = a + h;
    x2 = a + 2 * h;
    y0 = y(1);
    y1 = y(2);
    y2 = y(3);
    for k = 4 : n
        x = a+k*h;
        %Urmatoarele doua linii de cod se vor modifica!
        %Mai exact, se vor inlocui cu formulele date/deduse!
        ypr = y2 + h*( 23 * f(x2,y2)- 6 *f(x1,y1) +5 * f(x0,y0) )/12;
        ycor = ypr  +h*(5 * f(x,ypr) + 8 * f(x2,y2) - f(x1,y1) )/ 2;
        y(k) = ycor;
        y0 = y1;
        y1 = y2;
        y2 = ycor;
        x0 = x1;
        x1 = x2;
        x2 = x;
    endfor
endfunction



function y = Runge_Kutta(a, b, n, y0, f)
%   Functia care rezolva sistemul de ecuatii y' = f(x, y)
%prin metoda Runge - Kutta
%   Date de intrare:
%       - a, b -> intervalul de integrare;
%       - n    -> numarul de puncte;
%       - y0   -> conditia initiala;
%       - f    -> functia y'=f(x,y) .
%   Date de iesire:
%       - y    ->aproximarea solutiei în punctele xi.
    h = (b-a)/n;
    y = zeros(n+1,1);
    y(1) = y0;
    for k = 2 : n+1
        x = a + (k - 1)*h;
        K1 = h * f(x,  y(k-1) );
        K2 = h * f(x + h/2, y(k-1) + K1/2);
        K3 = h * f(x + h/2, y(k-1) + K2/2);
        K4 = h * f(x + h, y(k-1) + K3);
        y(k) = y(k-1) + (K1 + 2 * K2 + 2 * K3 + K4 )/6;
    endfor
endfunction