function cgmS

%  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,2
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];
        plot(1,1,'r*','MarkerSize',15,'LineWidth',3)
        hold on
    elseif ifun==2
        if ic==1
            x=[1, 0.8];
            exact=[5.535799358444e-01, -5.535799358444e-01];
        else
            x=[-1, -0.6];
            exact=[-4.187827176417e-01, 4.187827176417e-01];
        end
        xa=-1; xb=1; ya=-1; yb=1;
        xd=linspace(xa,xb,80);
        yd=linspace(ya,yb,80);
        [X,Y]=meshgrid(xd,yd);
        Z=(X-Y).^4+8*X.*Y-X+Y+3;
    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(['CGM-PR Step = 0','      Iterative Error = '],'FontSize',16,'FontWeight','bold')
    pause
    
    tol=0.000001;
    g=dF(ifun);
    d=-g;
    counter=0;
    while norm(x-exact)>tol & counter<40
        counter=counter+1;
        alpha=fminbnd(@F,0,0.5);
        x=x+alpha*d;
        g0=g;
        g=dF(ifun);
        beta=dot(g,g-g0)/dot(g0,g0);
        d=-g+beta*d;
        xpoints=[xpoints, x(1)];
        ypoints=[ypoints, x(2)];
        error=norm(alpha*d);
        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(['CGM-PR Step = ',num2str(counter),'      Iterative Error = ',num2str(error,'%2.1e')],'FontSize',16,'FontWeight','bold')
        pause
        
    end
    
    pause
    
end
fprintf('\n')

    function z=F(s)
        xs=x+s*d;
        if ifun==1
            z=100*(xs(1)^2-xs(2))^2+(xs(1)-1)^2;
        elseif ifun==2
            z=(xs(1)-xs(2))^4+8*xs(1)*xs(2)-xs(1)+xs(2)+3;
        end
    end

    function g=dF(ifun)
        if ifun==1
            g(1)=400*x(1)*(x(1)^2-x(2))+2*(x(1)-1);
            g(2)=-200*(x(1)^2-x(2));
        elseif ifun==2
            g(1)=4*(x(1)-x(2))^3+8*x(2)-1;
            g(2)=-4*(x(1)-x(2))^3+8*x(1)+1;
        end
    end

end