% Lichoběžníková metoda pro rovnici y'=-y^2

f=@(x,y)-y.^2;
x0=1;y0=1;
xn=5;

% Cyklus pro celý interval

h=0.5;
x=x0;
y=y0;
xx=x0;
yy=y0;
while x<xn
    P=[h/2, 1, y^2*h/2-y];
    r=roots(P);
    [m,I]=min(abs(r-y));
    y=r(I);
    x=x+h;
    xx=[xx,x];
    yy=[yy,y];
end
plot(xx,yy,xx,1./xx)

%
% Porovnání Impl. Eulerovy a Lichob. metody

f=@(x,y)-y.^2;
x0=1;y0=1;
xn=5;

% LM

h=0.5;
x=x0;
y=y0;
xx=x0;
yy=y0;
while x<xn
    P=[h/2, 1, y^2*h/2-y];
    r=roots(P);
    [m,I]=min(abs(r-y));
    y=r(I);
    x=x+h;
    xx=[xx,x];
    yy=[yy,y];
end

% IEM

x=x0;
y=y0;
XX=x0;
YY=y0;
while x<xn
    P=[h 1 -y];
    r=roots(P);
    [m,I]=min(abs(r-y));
    y=r(I);
    x=x+h;
    XX=[XX,x];
    YY=[YY,y];
end
plot(xx,yy,XX,YY,xx,1./xx,'--')