function  QIM_demo(f,I,nit)
% function  QIM_demo(f,I,nit)
% Quadrtic Interpolation Method demonstration
%  f  - function whose minimum is looked for (anonymouns function)
%  I  - inteval
% nit - number of iterations (default value is 5)

if nargin==2, nit=5; end
a=I(1); b=I(2); c=(a+b)/2;
fa=f(a); fb=f(b); fc=f(c);
x=linspace(a,b,201);
fx=f(x);
plot(x,fx,'b',[a b c],[fa fb fc],'r*');
t=title('Quadratic Interpolation Method Demonstration');

hold on
k=0;
err=0;
% nit=1;
xx=get(gca,'XLim');
yy=get(gca,'YLim');
fprintf('Iterations:\n')
while k<nit
    k=k+1;
    fab=(fb-fa)/(b-a);
    fbc=(fc-fb)/(c-b);
    fabc=(fbc-fab)/(c-a);
    P=[0 0 fa]+[0 fab*poly(a)]+fabc*poly([a,b]);
    Px=polyval(P,x);
    pause(1)
    plot(x,Px,'--')
    axis([xx yy])
    if fabc<=0, err=1; end
    d=(a+b-fab/fabc)/2;
    if d<a || d>b, err=1; end
    if err
	t1=text(xx(1),yy(1)-(yy(2)-yy(1))/10,'Parabola does not have a minimum in the interval.');
	set(t1,'Color',[1 0 0])
	break;
    else
        disp(d)
	fd=f(d);
	Pd=polyval(P,d);
	pause(1)
	plot(d,Pd,'ro')
	pause(1)
	plot(d,fd,'r*')
	if d>c
	    a=c; fa=fc;
	else
	    b=c;fb=fc;
	end
	c=d; fc=fd;
    end
end

end

