function xx=grad_demo(A,b,x0,ni)
% function xx=grad_demo(A,b,x0)
% Demonstration of gradient method

if nargin==3
    ni=20;
end
nx=101;
ny=101;
tol=1e-3;
x0=x0(:);
b=b(:);
s=A\b;
I=3*abs(x0-s);
I(I==0)=1;
xx=x0;
x=linspace(s(1)-I(1)/2,s(1)+I(1)/2,nx);
y=linspace(s(2)-I(2)/2,s(2)+I(2)/2,ny);
%[X,Y]=meshgrid(xx,yy);
f=zeros(ny,nx);
for ii=1:nx
    for jj=1:ny
        xy=[x(ii),y(jj)];
        f(jj,ii)=xy*A*xy'/2-xy*b(:);
    end
end
contour(x,y,f,20);
hold on;
fx=xx'*A*xx/2-xx'*b;
contour(x,y,f,[fx fx],'r--');
xn=xx;
xx=xn-1;

k=1;
while norm(xn-xx)>tol && (k<ni)
    xx=xn;
    s=b-A*xx;
    xn=xx+(s'*s)/(s'*A*s)*s;
    fx=xn'*A*xn/2-xn'*b;
    contour(x,y,f,[fx fx],'r--');
    plot([xx(1),xn(1)],[xx(2),xn(2)],'r');
    k=k+1;
    pause
end
hold off


end

