function fig9_5

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

clf
% get(gcf)
set(gcf,'Position', [24 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 (3 and 1 are in book)
iexample=3;

% set parameters
if iexample==1
    x=[1, 1];
    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',20,'FontWeight','bold')
    set(gca,'ytick',[-80 -60 -40 -20 0])
    nc=10;
    cx=linspace(cxL,cxR,nc);
    for i=1:nc
        cv(i)=0.5*(a*cx(i)^2 + 2*b*cx(i)*cy + d*cy*cy) - b1*cx(i) - b2*cy;
    end
elseif iexample==2
    x=[1, 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')
    ylabel('y-axis')
    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
elseif iexample==3
    x=[1, 2];
    a=3; b=-2; c=-2; d=4; b1=1; b2=-1;
    xL=-0.5; xR=2.5; yL=-0.5; yR=2.5;
    cxL=0; cxR=2; cy=1;
    axis([xL xR yL yR])
    xlabel('x-axis')
    ylabel('y-axis')
    cv=[-0.13  0.1 0.46 1.2 2.3 3.7 6.5 10];
    set(gca,'xtick',[-0.5 1.0 2.5])
    set(gca,'ytick',[-0.5 1.0 2.5])
elseif iexample==4
    x=[1, 1];
    a=2; b=1; c=1; d=3; b1=1; b2=-1;
    xL=0; xR=2; yL=-1; yR=1;
    cxL=0; cxR=2; cy=1;
    axis([xL xR yL yR])
    xlabel('x-axis')
    ylabel('y-axis')
    cv=[-0.68 -0.585 -0.3 0.4 1.3 2.3 3.5 4.8];
end
tol=0.0001;

% generate surface data for plotting
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;
% surf(X,Y,Z)
% pause

% 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;
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\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.6)
    plot(xpoints,ypoints,'r-','LineWidth',2.5)
    plot(xpoints,ypoints,'o','Color',[110/256 0 0],'MarkerSize',15,'LineWidth',2.5)
    error=norm(diff);
    set(gca,'FontSize',20,'FontWeight','bold')
    hold off
end

%0.5*(a*x(1)^2 + 2*b*x(1)*x(2) + d*x(2)*x(2)) - b1*x(1) - b2*x(2)

xpoints(2)
ypoints(2)

A=[[a b];[c d]];
xe=A\[b1 b2]';
length(xpoints)
error=norm([x(1) x(2)]'-xe,inf)
eigenvalues=eig(A)
cond(A,inf)
cond(A,2)


% if iexample==3
%     exportgraphics(gcf,'/Users/mark/Desktop/sdm01.eps')
% elseif  iexample==1
%     exportgraphics(gcf,'/Users/mark/Desktop/sdm02.eps')
% end





