function fig9_17

% i_fun = number of function evaluations
global i_fun

%  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

i_fun=0;

x=[3/4   -3/4]';
exact=[1, 1]';
tol=1e-6;
xpoints=x(1);
ypoints=x(2);
f=F(x);
ff=norm(f,2)^2;
D=zeros(2,2);
mu=10;
alpha=2;
beta=2;
counter=1;
err=10*tol;
%while err>tol %& counter<200
while norm(x-exact,inf)>tol
    counter=counter+1;
    xx=x;
    ff0=ff;

    J=Jac(x);
    JJ=J'*J;
    gk=2*J'*f;
    D(1,1)=2*JJ(1,1);  D(2,2)=2*JJ(2,2);
    S=2*JJ+mu*D;

    % check on descent direction
    %eig(S)

    dx=-inv(S)*gk;
    x=xx+dx;
    f=F(x);
    ff=norm(f,2)^2;

    if ff<ff0
        mu=mu/beta;
        xpoints=[xpoints, x(1)];
        ypoints=[ypoints, x(2)];
        err=norm(inv(S)*gk./x,inf);

        % check on jacobian approximation
        fx=20*x(1); fy=-10; gx=1; gy=0;
        fxx=20; fxy=0; fyy=0; gxx=0; gxy=0; gyy=0;
        [(f(1)*fxx+f(2)*gxx)/(fx^2+gx^2) (f(1)*fxy+f(2)*gxy)/(fx*fy+gx*gy) (f(1)*fyy+f(2)*gyy)/(fy^2+gy^2)]

    else
        mu=alpha*mu; x=xx;  ff=ff0;

    end

end

i_fun
% check on jacobian approximation

clf
hold on
% get(gcf)
set(gcf,'Position', [12 950 421 395])
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)

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;
plot(exact(1),exact(2),'r*','MarkerSize',14,'LineWidth',2)

box on
contour(X,Y,Z,40,'b','LineWidth',1)
axis([xa-0.004 xb ya yb])
axis square
plot(xpoints,ypoints,'r-','LineWidth',2)
plot(xpoints,ypoints,'o','Color',[100/256 0 0],'MarkerSize',10,'LineWidth',2)

set(gca,'xtick',[-1 -0.5 0 0.5 1])
set(gca,'ytick',[-1 -0.5 0 0.5 1])

xlabel('x-axis')
ylabel('y-axis')
set(gca,'FontSize',16,'FontWeight','bold')
hold off

%exportgraphics(gcf,'/Users/mark/Desktop/LMM.eps')

% F=dot(f,f)
function f=F(x)
global i_fun
i_fun=i_fun+1;
xx=x(1); yy=x(2);
f1=10*(xx^2-yy);
f2=xx-1;
f=[f1;f2];

% jacobian
function g=Jac(x)
global i_fun
i_fun=i_fun+1;
xx=x(1); yy=x(2);
dfx=20*xx;  dfy=-10;
dgx=1; dgy=0;
g=[dfx,dfy;dgx,dgy];




















