function [xx,yy]=lm_demo(f,x0,xn,y0,h,u)

x=x0;
y=y0;
xx=x0;
yy=y0;
ff=strrep(f,'*','.*');
ff=strrep(ff,'/','./');
ff=strrep(ff,'^','.^');
fy=diff(str2sym(f),'y');
while x<xn
    fi=eval(f);
    yi=y;
    xi=x;
    x=x+h;
    z0=y+h*fi;
    y=z0;
    z1=z0-(z0-yi-0.5*h*(fi+eval(f)))/(1-0.5*h*eval(fy));
    while abs(z1-z0)>0.001
        y=z1;
        z0=z1;
        z1=z0-(z0-yi-0.5*h*(fi+eval(f)))/(1-0.5*h*eval(fy));
    end
    y=z1;
    xx=[xx,x];
    yy=[yy,y];
end
miny=min(yy);
maxy=max(yy);
% figure
hold on
plot(x0,y0,'r*');
axis([x0,xn,miny,maxy])
t1=title(['Lichobeznikova metoda, h=',num2str(h),', rovnice y=',f,]);
set(t1,'fontsize',16);
for k=2:length(xx)
    line([xx(k-1),xx(k)],[yy(k-1),yy(k)]);
    plot(xx(k),yy(k),'r*')
    pause(0.1)
end
if nargin==6
    pomx=linspace(x0,xn);
    x=pomx;
    pomy=eval(u);
    plot(pomx,pomy,'--g')
end