clear all; close all;
f = @(x) (x(1)+sqrt(x(2))).^2+10.*cos(x(1).^2+x(2)).^2;
xlim=[-2,0 ;0,5];
h = 10^-8;
h1 = [h 0]; h2 = [0 h];
path=[];
a = 30;
%for k = 1:10
%    a=a/k;
%    a
x0  = [-1 1.5];x=x0;
path=[path ;x0];
n = 0;


fx = f(x);f0=fx+10;
fprintf('%d %20.15f %20.15f %20g\n',n,x,fx)
n=n+1;
while abs(fx-f0)>10^-16
    f0 = fx;
    % grad f
    gf = @(x) [f(x+i*h1) f(x+i*h2)]/h;
    gg = imag(gf(x));
    % Hessian
    g1 = imag(gf(x+h1)-gf(x-h1))/(2*h);
    g2 = imag(gf(x+h2)-gf(x-h2))/(2*h);
    H= [g1 ; g2 ];
    H = H+a.*eye(size(H,1),size(H,2));
    p  = linsolve(H,-gg');
    x  = x+p';    
    path=[path;x];
    fx = f(x);
    fprintf('%d %20.15f %20.15f %20g\n',n,x,fx)
    n = n+1;
    if n>100 
        break;
    end
end
%end
x=-2:0.1:0;
y=0:0.1:5;
[X,Y]=meshgrid(x,y);
for i=1:size(X,1)
        for j=1:size(X,2)
            x=[X(i,j) Y(i,j)];
            Z(i,j)=f(x);
        end
end
contour(X,Y,Z,10) ;colorbar;hold on;
x=[path(:,1) path(:,2)];
plot3(path(:,1),path(:,2),zeros(size(path,1)),'mo-');
ende=plot3(path(size(path,1),1),path(size(path,1),2),0,'bo','MarkerSize',10);
legend(ende,'final');

