Оптимизация функции двух переменных методом покоординатного спуска

clear all; clc;
format long;

g = 0.5; % постоянная шага
d = 0.01; % дельта
% Начальная точка
x1 = 0; 
x2 = 0; 
k = 1;  % Счетчик шагов
kmax = 1000; % Предельное число шагов, 
% задается для предотвращения зацикливания
% Массивы для хранения промежуточных координат
x1trace = [x1]; 
x2trace = [x2]; 
i = 2; 
while k < kmax 
        % Спуск по первой координате
        gr1 = 6*x1 - 2*x2 + 1;  % производная по х1
        x1 = x1 + g*gr1;
        % Сохранение координат
        x1trace(i) = x1; 
        x2trace(i) = x2; 
        i = i + 1; 
        % Спуск по второй координате
        gr2 = 8*x2 - 2*x1 - 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 = -100:0.1:100; 
y = -100:0.1:100; 
[X, Y] = meshgrid(x, y); 
Z=3*X.^2-2*X.*Y+4.*Y.^2+X-Y;        % исходная функция
[C, h] = contour(X, Y, Z); 
clabel(C, h); 
hold on;
plot(x1trace, x2trace, '-'); 
% Вывод начальной точки на график
text(x1trace(1) + 0.2, x2trace(1) + 0.5, 'M0'); 
% Вывод решения на график
text(x1 + 2, x2, ...
strvcat(['x1 = ' (num2str(x1))], ...
        ['x2 = ' (num2str(x2))], ...
        ['k = '  (num2str(k))])); 
hold off;

Проблема в том, что программа упирается в ограничитель итераций, даже если поставить ограничение 10000. Производную функции я искал через diff и потом проверил, в этом ошибки быть не должно. Код программы я брал готовый, просто заменил функцию и ее производные на свои. Со старой функцией программе хватало 50 итераций.

С начальными условиями, точностью и шагом экспериментировал. Результат тот же.