function sdm11

%  SDM for 2x2 system
%        [a b][x] = [b1]
%        [c d][y] = [b2]


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)

% specify example
iexample=2

% 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])
    xlabel('x-axis','FontSize',24,'FontWeight','bold')
    ylabel('y-axis','FontSize',24,'FontWeight','bold')
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])
    xlabel('x-axis','FontSize',24,'FontWeight','bold')
end
A=[[a b];[c d]]
cond_A=cond(A,2)
tol=0.01;
x=[1, 1];

% plot surface for quadratic form
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)

pause

clf
set(gcf,'Position', [1 868 482 477])
% determine 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)];
diff=[1,1];
counter=1;
hold on
box on
axis square
contour(X,Y,Z,cv,'b','LineWidth',1.1)
plot(xpoints,ypoints,'ko','MarkerSize',10,'LineWidth',2)
xlabel('x-axis')
ylabel('y-axis')
set(gca,'FontSize',24,'FontWeight','bold')

pause

while norm(diff)>tol & (counter<30)
    counter=counter+1;
    q=[a*r(1)+b*r(2),c*r(1)+d*r(2)];
    alpha=dot(r,r)/dot(r,q);
    diff=alpha*r;
    x=x+diff;
    r=r-alpha*q;
    
    fprintf('\n  %i  Computed Solution =  %e  %e \n',counter,x)
    
    xpoints=[xpoints, x(1)];
    ypoints=[ypoints, x(2)];
    hold on
    box on
    axis square
    contour(X,Y,Z,cv,'b','LineWidth',1.1)
    plot(xpoints,ypoints,'r-','LineWidth',2)
    plot(xpoints,ypoints,'ko','MarkerSize',10,'LineWidth',2)
    error=norm(diff);
    set(gca,'FontSize',24,'FontWeight','bold')
    hold off
    
    pause
    
end
fprintf('\n',counter,x)