function cgm1

%  CGM for 2x2 system (c=b, a>0, and ad-b^2>0)
%        [a b][x1] = [b1]
%        [c d][x2] = [b2]


% specify example
iexample=2
tol=0.01;

% set parameters
if iexample==1
    a=2; b=2; c=b; d=3; b1=1; b2=-1;
    xL=0; xR=6; yL=-5; yR=1;
    cxL=3; cxR=10; cy=-2;
    axis([0 4 -3 1])
elseif iexample==2
    a=5; b=4.99; c=b; d=5; b1=1; b2=-1;
    xL=-2; xR=80; yL=-80; yR=5;
    cxL=0; cxR=100; cy=-100;
    axis([xL xR yL yR])
end

A=[[a b];[c d]]
cond_A=cond(A,2)
x=[1, 1];

% plot surface for quadratic form
clf
% get(gcf)
set(gcf,'Position', [1 868 482 477])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
xd=linspace(xL,xR,60);
yd=linspace(yL,yR,60);
[X,Y]=meshgrid(xd,yd);
Z = 0.5*(a*X.^2 + 2*b*X.*Y + d*Y.*Y) - b1*X - b2*Y;
surfc(X,Y,Z)
ylabel('y-axis')
xlabel('x-axis')
set(gca,'FontSize',18,'FontWeight','bold')

pause

clf
set(gcf,'Position', [1 868 482 477])
% specify elevations to draw contours
cx=linspace(cxL,cxR,18);
for i=1:18
    cv(i)=0.5*(a*cx(i)^2 + 2*b*cx(i)*cy + d*cy*cy) - b1*cx(i) - b2*cy;
end

% start iteration
xpoints=x(1);
ypoints=x(2);
r=[b1-a*x(1)-b*x(2),b2-c*x(1)-d*x(2)];
dd=r;
diff=[1,1];
counter=0;
error=norm(diff);

contour(X,Y,Z,cv,'b','LineWidth',1.5)
hold on
box on
axis([xL xR yL yR])
axis square
plot(xpoints,ypoints,'ko','MarkerSize',12,'LineWidth',2)
ylabel('y-axis')
xlabel('x-axis')
set(gca,'FontSize',18,'FontWeight','bold')
title(['Step = 0','      Iterative Error = '],'FontSize',16,'FontWeight','bold')

pause

while norm(diff)>tol
    counter=counter+1;
    if counter==1
        beta=0;
    else
        beta=dot(r,r)/dot(r0,r0);
    end
    dd=r+beta*dd;
    q=[a*dd(1)+b*dd(2),c*dd(1)+d*dd(2)];
    alpha=dot(r,r)/dot(dd,q);
    diff=alpha*dd;
    x=x+diff;
    r0=r;
    r=r-alpha*q;
    
    fprintf('\n  %i  Computed Solution =  %e  %e \n',counter,x);
    
    xpoints=[xpoints, x(1)];
    ypoints=[ypoints, x(2)];
    shading interp
    hold on
    box on
    contour(X,Y,Z,cv,'b','LineWidth',1.5)
    plot(xpoints,ypoints,'r-','LineWidth',2)
    plot(xpoints,ypoints,'ko','MarkerSize',10,'LineWidth',2)
    error=norm(diff);
    set(gca,'FontSize',18,'FontWeight','bold')
    title(['Step = ',num2str(counter),'      Iterative Error = ',num2str(error)],'FontSize',16,'FontWeight','bold')
    
    pause
    
end
fprintf('\n')