function [xx,yy] = PK_multistep(f,I,h,y0)
% function [xx,yy] = PK_multistep(f,I,h,y0)
% vicekrokova metoda prediktor-koerktor
% prediktor: AB4
% korektor: Simpsonova metoda
%   
presnost=1e-12;
a=[1 0 0 0];
b=[55/24 -59/24 37/24 -9/24];
la=length(a);           % pocet kroku
x0=I(1);                % pocatecni uzel
xn=I(2);                % koncovy uzel
x1=x0+h*(0:la-1);       % body pro poc. hodnoty
[xx,yy]=ode45(f,x1,y0); % poc. hodnoty
if la==2                % metoda je dvoukrokova
    xx=[xx(1) xx(end)];
    yy=[yy(1) yy(end)];
else
    xx=xx(:)';          % vysledky dame do radku
    yy=yy(:)';
end
yd=f(xx,yy);            % priblizne hodnoty derivaci
nx=length(xx);          % pocet urcenych hodnot
x=xx(nx);               % posledni spocitany uzel
while x<xn              % hlavni cyklus vypoctu
    I=nx:-1:nx-la+1;    % indexy pouzitych uzlu pozpatku
    y=sum(a.*yy(I))+h*sum(b.*yd(I));    % vicekrokova formule - prediktor AB4
    x=x+h;              % dalsi uzel
    d=f(x,y);           % derivace v dalsim uzlu
    y1=yy(nx-1)+h*(yd(nx-1)+4*yd(nx)+d)/3; % korektor - Simpsonova metoda
    while abs(y1-y)>presnost
        y=y1;
        d=f(x,y);
        y1=yy(nx-1)+h*(yd(nx-1)+4*yd(nx)+d)/3; % korektor - Simpsonova metoda
    end
    y=y1;
    d=f(x,y);
    yy=[yy,y];          % doplneni hodnot
    xx=[xx,x];
    yd=[yd,d];
    nx=nx+1;            % zvysime posledni index
end
end

