function sdmS

%  Inputs:
%        x = starting vector
%        xa, xb = x-interval used in contour plot
%        ya, yb = y-interval used in contour plot
%        tol = tolerance for stopping iteration

%  Required (end of file)
%        F    is function F(x,y) to be minimized
%        dF  is gradient of F(x,y)

%  Specify which example to run
ifun=1


figure(1)
clf
hold on
% get(gcf)
set(gcf,'Position', [4 431 560 420])
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)

for ic=1:2
    
    if ifun==1
        if ic==2
            return
        end
        x=[-1.0, -0.5];
        %x=[-1.2, 1.8];
        xa=-1.1; xb=1.1; ya=-1.1; yb=1.1;
        xd=linspace(xa,xb,80);
        yd=linspace(ya,yb,80);
        [X,Y]=meshgrid(xd,yd);
        Z=100*(X.^2-Y).^2+(X-1).^2;
        exact=[1, 1];
        tol=1e-7;
        plot(1,1,'r*','MarkerSize',15,'LineWidth',3)
        hold on
    elseif ifun==2
        if ic==1
            x=[1, 0.8];
            exact=[.55357993584438379685, -.55357993584438379685];
        else
            x=[-1, -0.6];
            exact=[-4.187827176417e-01, 4.187827176417e-01];
        end
        xa=-1; xb=1; ya=-1; yb=1;
        nc=80;
        xd=linspace(xa,xb,nc);
        yd=linspace(ya,yb,nc);
        [X,Y]=meshgrid(xd,yd);
        Z=(X-Y).^4+8*X.*Y-X+Y+3;
        tol=1e-4;
    elseif ifun==3
        x=[1, 1];
        xa=-1; xb=1; ya=-1; yb=1;
        xd=linspace(xa,xb,80);
        yd=linspace(ya,yb,80);
        [X,Y]=meshgrid(xd,yd);
        Z=10*X.^2+Y.^2;
        exact=[0, 0];
        tol=1e-7;
    end
    
    xpoints=x(1);
    ypoints=x(2);
    box on
    contour(X,Y,Z,40,'b','LineWidth',1)
    axis([xa-0.004 xb ya yb])
    axis square
    xlabel('x-axis')
    ylabel('y-axis')
    plot(xpoints,ypoints,'ko','MarkerSize',12,'LineWidth',3)
    set(gca,'xtick',[-1 -0.5 0 0.5 1])
    set(gca,'ytick',[-1 -0.5 0 0.5 1])
    set(gca,'FontSize',18,'FontWeight','bold')
    title(['SDM Step = 0','      Iterative Error = '],'FontSize',16,'FontWeight','bold')
    pause
    
    [dfx, dfy]=feval('dF',x(1),x(2),ifun);
    df1=[dfx, dfy];
    d1=-df1;
    q=1;
    xx=10;
    counter=0;
    
    while norm(xx-exact)>tol & counter<50
        counter=counter+1;
        options = optimset('TolX',tol);
        s=fminbnd(@F,0,1,options,x(1),x(2),d1(1),d1(2),ifun);
        q=s*d1;
        xx=x+q;
        [dfx, dfy]=feval('dF',xx(1),xx(2),ifun);
        df1=[dfx, dfy];
        d1=-df1;
        x=xx;
        xpoints=[xpoints, x(1)];
        ypoints=[ypoints, x(2)];
        error=norm(q);
        plot(xpoints,ypoints,'r-','LineWidth',1.5)
        plot(xpoints,ypoints,'ko','MarkerSize',10,'LineWidth',3)
        
        fprintf('\n  %i  Computed Solution =  %e  %e \n',counter,x)
        title(['SDM Step = ',num2str(counter),'      Iterative Error = ',num2str(error,'%2.1e')],'FontSize',16,'FontWeight','bold')
        pause
        
    end
    
    pause
    
end




function z=F(s,xx,yy,dx,dy,ifun)
x=xx+s*dx; y=yy+s*dy;
if ifun==1
    z=100*(x^2-y)^2+(x-1)^2;
elseif ifun==2
    z=(x-y)^4+8*x*y-x+y+3;
elseif ifun==3
    z=10*x^2+y^2;
end


function [dfx, dfy]=dF(x,y,ifun)
if ifun==1
    dfx=400*x*(x^2-y)+2*(x-1);
    dfy=-200*(x^2-y);
elseif ifun==2
    dfx=4*(x-y)^3+8*y-1;
    dfy=-4*(x-y)^3+8*x+1;
elseif ifun==3
    dfx=20*x;
    dfy=2*y;
end