clear
close all
R = 0.125;
epsilon = 0.001;
x0=[1 4];
[x,y]=meshgrid(-10:.2:10);
z = (x + 2).^2 + (y + 2).^2 + ((y - x).^2)/R;
hold on
contour(x,y,z);
[dx,dy]=gradient(z);
quiver(x,y,dx,dy);
xs=x0;
itsk = 100;
step=0.05;
   for i=1:itsk
        x0 = xs - step * gradient1(xs,R);
        plot([xs(1),x0(1)],[xs(2),x0(2)],'r',[xs(1),x0(1)],[xs(2),x0(2)],'ro')
        d = abs(xs - x0);
        xs=x0;
       % R = R/5;
        if(d < epsilon)
        break;
        end
   end
disp(x0);
function gr = gradient1(x,R)
gr(1) = (2*(R*(x(1) + 2) + x(1) - x(2)))/R;
gr(2) = (2*(R*(x(2) + 2) - x(1) + x(2)))/R;
end