function [xx,yy]=lm(f,x0,xn,y0,h)
% function [xx,yy]=lm(f,x0,xn,y0,h)
% lichobeznikova  metoda pro rovnici y'=f(x,y)
% f  - anonymni funkce promennych x,y
% x0 - zacatek intervalu
% xn - konec intervalu
% y0 - pocatecni hodnota
% h  - krok
% nelinearni rovnice se resi Newtonovou metodou

presnost=0.000001;
x=x0;
y=y0;
xx=x0;
yy=y0;
fs_y=diff(sym(f),'y');		                % derivace f podle y symbolicky
fy=matlabFunction(fs_y,'Vars',{'x','y'});   % derivace f podle y
while x<xn
    fi=f(x,y);      % hodnota funkce f v predchoyim bode
    z0=y+h*f(x,y);
    yi=y;
    xi=x;
    x=x+h;
    z0=y+h*fi;      % pocatecni aproximace reseni nelin. rovnice je hodnota ziskana s pomoci EEM
    y=z0;
    z1=z0-(z0-yi-0.5*h*(fi+f(x,y)))/(1-0.5*h*fy(x,y));    % prvni aproximace rseni nelin. rovnice
    while abs(z1-z0)>0.001                                  % Newtonova metoda
        y=z1;
        z0=z1;
        z1=z0-(z0-yi-0.5*h*(fi+f(x,y)))/(1-0.5*h*fy(x,y));
    end
    y=z1;
    xx=[xx,x];
    yy=[yy,y];
end
