function x1=NS_demo(A,b,x0,ni)
% function NS_demo(A,b,x0,ni)
% demonstrace metody nejvetsiho spadu pro system rovnic Ax=b
% x0 - pocatecni iterace
% ni - pocet iteraci (implicitne 10)
%

if nargin==3
    ni=10;
end
[m,n]=size(A);
[V,D]=eig(A);
D=diag(D);
if m~=2 || n~=2 || any(any(A'~=A)) || any(D<=0)
    error('Matice musi byt pozitivne definitni 2x2!');
end
b=b(:);x0=x0(:); % prevod na sloupce
[m,n]=size(b);
if m~=2 || n~=1
    error('Spatny rozmer vektoru b!');
end
[m,n]=size(x0);
if m~=2 || n~=1
    error('Spatny rozmer vektoru pocatecniiterace x0!');
end

Q=@(x,y)[x y]*A*[x; y]/2-[x y]*b;

s=inv(A)*b;
Qs=Q(s(1),s(2));
Q0=Q(x0(1),x0(2));
Qs0=Qs+0.0004*(Q0-Qs);
d=sqrt(2*(Q0-Qs)/min(D));
LList=[Qs0,Q0];
warning ('off','all');
Qcont=fcontour(Q,[s(1)-1*d,s(1)+1*d,s(2)-1*d,s(2)+1*d]);
%Qcont=fcontour(Q);
set(Qcont,'LevelList',LList,'LineColor','green','LineStyle','--');
hold on
plot(x0(1),x0(2),'r*');
r=b-A*x0;
for k=1:ni
%    title('Neco stiskni.')
    pause(3)
    t=r'*r/(r'*A*r);
    x1=x0+t*r;
    Q1=Q(x1(1),x1(2));
    LList=[LList,Q1];
    set(Qcont,'LevelList',LList);
    plot(x1(1),x1(2),'r*')
    plot([x0(1),x1(1)],[x0(2),x1(2)],'b--');
    r=r-t*A*r;
    x0=x1;
end
title('HOTOVO!')
hold off
