%  quadratic case – conjugate gradient
clear
itmax = 2;
X= [1 1; 0 1];Y=[-1; 1];

for i=1:41
   xx(i)=-4+(i-1)*0.1;
   for j=1:41
      yy(j)=-1+(j-1)*0.1;
      BB=[xx(i);yy(j)];
      Z(i,j) = (Y-X*BB)'*(Y-X*BB);
   end
end
figure(1); clf;
contour(xx,yy,Z,20)
xlabel('B_1')
ylabel('B_2')
hold on
clear xx; clear yy   
xtickformat('%0.1f')
ytickformat('%0.1f') 


B=[0;0];xx(1)=0;yy(1)=0;g=0;
for k=1:itmax
   g=-X'*(Y-X*B);ng=g'*g;
   if k==1
      w=-g;
   else
      w=-g+ng/ng1*w1;
   end
   ro=-(g'*w)/(w'*X'*X*w)
   B=B+ro*w
   g1=g;w1=w;ng1=ng;
yy(k+1)=B(2);xx(k+1)=B(1);
s(k)=0.5*(Y-X*B)'*(Y-X*B)
end
plot(xx,yy,'lineWidth',2)

figure(2); clf;
plot(s, 'o-', 'linewidth', 1.0)
xlabel('iteration');
ylabel('S=(Y-XB)^T(Y-XB)')
ytickformat('%0.2f')
xticks( [1 2] )
grid on


% fig=figure(2);
% print(fig,fullfile(pwd, 'Figure6_12.tif'),'-dtiff','-r600')
% 
% fig=figure(1)
% print(fig,fullfile(pwd, 'Figure6_13.tif'),'-dtiff','-r600')