Определить минимум функции по методу Гаусса-Зейделя

% Метод Гаусса-Зейделя.
clear all;
clc;
% Значения коэффициентов.
A1 = 3;
A2 = 2;
a1 = 4;
a2 = 2;
 b = 1;
c1 = 4;
c2 = 2;
D = 2;
g = 0.2; % Постоянная шага.
d = 0.01; % дельта.
% Начальная точка.
x1 =5 ;
x2 =8;
k = 1;  % Счетчик шагов
kmax = 100; % Предельное число шагов,
% задается для предотвращения зацикливания.
% Массивы для хранения промежуточных координат.
x1trace = [x1];
x2trace = [x2];
i = 2;
while k < kmax
% Спуск по первой координате.
gr1 = -A1*(2*a1-2*x1)/b.^2*(((a1-x1).^2/b.^2)+((c1-x2).^2/D.^2)+1).^2-A2*(2*a2-2*x1)/b.^2*(((a2-x1).^2/b.^2)+((c2-x2).^2/D.^2)+1).^2;
x1 = x1 + g*gr1;
% Сохранение координат.
x1trace(i) = x1;
x2trace(i) = x2;
i = i + 1;
% Спуск по второй координате.
gr2 = -A1*(2*c1-2*x2)/D.^2*(((a1-x1).^2/b.^2)+((c1-x2).^2/D.^2)+1).^2-A2*(2*c2-2*x2)/D.^2*(((a2-x1).^2/b.^2)+((c2-x2).^2/D.^2)+1).^2;
x2 = x2 + g*gr2;
% Сохранение координат.
x1trace(i) = x1;
x2trace(i) = x2;
i = i + 1;
% Проверка условия останова.
if sqrt(gr1^2 + gr2^2) <= d;
break; % Выход из цикла в случае выполнения условия.
end
k = k + 1;
end
% Построение графика.
x = -10:0.1:10;
y = -10:0.1:10;
[X, Y] = meshgrid(x, y);
Z=100-A1./(1+((X-a1)./b).^2+((Y-c1)./D).^2)-A2./(1+((X-a2)./b).^2+((Y-c2)./D).^2);
[C, h] = contour(X, Y, Z);
clabel(C, h);
% Отображение меток на линиях уровня.
hold on;
plot(x1trace, x2trace, '-');
grid on;
% Вывод начальной точки на график.
text(x1trace(1) + 0.2, x2trace(1) + 0.5, 'M0');
% Вывод решения на график.
text(x1 + 2, x2, ...
strvcat(['x1 = ' (num2str(x1))], ...
        ['x2 = ' (num2str(x2))], ...
        ['k = '  (num2str(k))]));

inf Infinity.       
    inf returns the IEEE arithmetic representation for positive
    infinity.  Infinity is also produced by operations like dividing by
    zero, eg. 1.0/0.0, or from overflow, eg. exp(1000).

syms x y A1 A2 a1 a2 b c1 c2 D
z = 100-A1/(1+((x-a1)/b)^2+((y-c1)/D)^2)-A2/(1+((x-a2)/b)^2+((y-c2)/D)^2)
diff(z, x)
diff(z, y)

% Метод Гаусса-Зейделя.
% Значения коэффициентов.
A1 = 3;
A2 = 2;
a1 = 4;
a2 = 2;
 b = 1;
c1 = 4;
c2 = 2;
D = 2;
g = 0.2; % Постоянная шага.
d = 0.01; % дельта.
% Начальная точка.
x1 =5 ;
x2 =8;
k = 1;  % Счетчик шагов
kmax = 100; % Предельное число шагов,
% задается для предотвращения зацикливания.
% Массивы для хранения промежуточных координат.
x1trace = [x1];
x2trace = [x2];
i = 2;
while k < kmax
% Спуск по первой координате.
% gr1 = -A1*(2*a1-2*x1)/b.^2*(((a1-x1).^2/b.^2)+((c1-x2).^2/D.^2)+1).^2-A2*(2*a2-2*x1)/b.^2*(((a2-x1).^2/b.^2)+((c2-x2).^2/D.^2)+1).^2;
gr1 = - (A1*(2*a1 - 2*x1)) / (b^2*((a1 - x1).^2/b^2 + (c1 - x2).^2/D^2 + 1).^2) - (A2*(2*a2 - 2*x1))/(b^2*((a2 - x1)^2/b^2 + (c2 - x2)^2/D^2 + 1)^2);
x1 = x1 - g*gr1;
% Сохранение координат.
x1trace(i) = x1;
x2trace(i) = x2;
i = i + 1;
% Спуск по второй координате.
% gr2 = -A1*(2*c1-2*x2)/D.^2*(((a1-x1).^2/b.^2)+((c1-x2).^2/D.^2)+1).^2-A2*(2*c2-2*x2)/D.^2*(((a2-x1).^2/b.^2)+((c2-x2).^2/D.^2)+1).^2;
gr2 = - (A1*(2*c1 - 2*x2))/(D^2*((a1 - x1)^2/b^2 + (c1 - x2)^2/D^2 + 1)^2) - (A2*(2*c2 - 2*x2))/(D^2*((a2 - x1)^2/b^2 + (c2 - x2)^2/D^2 + 1)^2);
x2 = x2 - g*gr2;
% Сохранение координат.
x1trace(i) = x1;
x2trace(i) = x2;
i = i + 1;
% Проверка условия останова.
if sqrt(gr1^2 + gr2^2) <= d;
break; % Выход из цикла в случае выполнения условия.
end
k = k + 1;
end
% Построение графика.
x = -10:0.1:10;
y = -10:0.1:10;
[X, Y] = meshgrid(x, y);
Z=100-A1./(1+((X-a1)./b).^2+((Y-c1)./D).^2)-A2./(1+((X-a2)./b).^2+((Y-c2)./D).^2);
[C, h] = contour(X, Y, Z);
clabel(C, h);
% Отображение меток на линиях уровня.
hold on;
plot(x1trace, x2trace, '-');
grid on;
% Вывод начальной точки на график.
text(x1trace(1) + 0.2, x2trace(1) + 0.5, 'M0');
% Вывод решения на график.
text(x1 + 2, x2, ...
strvcat(['x1 = ' (num2str(x1))], ...
        ['x2 = ' (num2str(x2))], ...
        ['k = '  (num2str(k))]));