function [u,k,konvergHist] = metodaNejvetsihoSpadu(A,b,epsil)

n = length(b);
u = zeros(n,1);
r = b-A*u;
normr0 = norm(r);
normr = normr0;
k = 0;
konvergHist(k+1) = normr;
Uk = u;
while normr/normr0 > epsil
    alfa = normr^2/(r'*A*r);
    u = u + alfa*r;
    Uk = [Uk,u];
    r = b - A*u;
    normr = norm(r);
    k = k+1;
    konvergHist(k+1) = normr;
end

if n==2
    upresne = A\b;
    u1 = linspace(min(Uk(1,:)),max(Uk(1,:)),100);
    u2 = linspace(min(Uk(2,:)),max(Uk(2,:)),100);
    [U1,U2] = meshgrid(u1,u2);
    Phi = 0.5*(A(1,1)*U1.^2+A(1,2)*U1.*U2+A(2,1)*U2.*U1+A(2,2)*U2.^2) - b(1)*U1-b(2)*U2;
    contour(U1,U2,Phi);
    axis equal
    
    hold on
    plot(upresne(1),upresne(2),'r*')
    hold on; plot(Uk(1,:),Uk(2,:),'b')
end