Оптимизация методом Ньютона (нахождение точки минимума). Оптимизация кода

%Метод Ньютона
clear; clc; clf;
format;
 
[X1, X2] = meshgrid(-2.5 : .1: 2.5, -2.5 : .1 : 2.5); % Сетка для линий уровня
 
% Исходные данные
f = '(2*X1^2-X1*X2+3*X2^2-X1+2*X2)'; % задаём функцию
FX = 2*X1.^2 - X1.*X2 + 3*X2.^2 - X1 + 2*X2; % функция для линий уровня
XK = [2, 1]; % Задаём X0
 
% Находим производные
dfX1 = char(diff(str2sym(f), ('X1'))); % частная производная по X1
dfX2 = char(diff(str2sym(f), ('X2'))); % частная производная по X2
d2fX1 = char(diff(str2sym(f), ('X1'), 2)); % вторая производная по X1
d2fX1X2 = char(diff(str2sym(dfX1), ('X2'))); % вторая производная по X1,X2
d2fX2 = char(diff(str2sym(f), ('X2'), 2)); % вторая производная по X2
F = inline(f); % преобразование
FX1 = inline(dfX1);
FX2 = inline(dfX2);
F2X1 = inline(d2fX1);
F2X12 = inline(d2fX1X2);
F2X2 = inline(d2fX2);
 
% Подготавливаем данные для запуска цикла приближений
x1 = {0, XK(1)}; % Формируем массив X1 координат, заносим Х0
x2 = {0, XK(2)}; % Формируем массив X2 координат, заносим Х0
ZK = F((XK(1)),(XK(2))); %Находим значение функции с текущим XK
ZK1 = ZK + 1; % Искусственно завышаем начальное значение функции для запуска цикла
E = 0.001; % задаем точность
N = 0; % Инициализация счетчика
HX = [F2X1(0) F2X12(0); F2X12(0) F2X2(0)]; % матрица Гессе (Нули, т.к. нет переменных)
IH = inv(HX); % Обратная матрица
 
% Запускаем цикл приближений
while abs(ZK1 - ZK) > E
    gradfXK = [FX1(XK(1),XK(2)); FX2(XK(1),XK(2))]; % формируем градиент
    ZK = F((XK(1)),(XK(2))); %Находим значение функции с текущим XK
    XK1 = XK; % Заносим предыдущее значение Х (Х0 на первой итерации)
    XK = XK1 - IH*gradfXK; % Находим X*
    ZK1 = F((XK(1)),(XK(2))); % Находим значение функции с найденной точкой
    if ZK1 < ZK, XK1 = XK; end % Проверяем, меньше ли значение функции
    if ZK1 ~= ZK, N = N + 1; end % Если приближение свершилось, увеличиваем счетчик
    if ZK1 == ZK continue;end % Если нет, то координаты точки в массив не добавляем
    x1(end+1,:) = {N, XK(1)}; % добавляем в массив координаты
    x2(end+1,:) = {N, XK(2)}; % найденной на данной итерации точки
end
 
% Вывод результатов
A = ['X* ='];
disp(A);
K = [XK(1); XK(2)];
disp(K)
M = ['Значение найдено на ',num2str(N),'-ой итерации.'];
disp(M)
 
% Линии уровня
x1 = cell2mat(x1(:,2)); % Извдекаем координаты
x2 = cell2mat(x2(:,2)); % точек из массивов
[CMatr, h] = contour(X1,X2,FX, 10);
clabel(CMatr, h);
grid on;
hold on;
plot(x1,x2,'r-o');
 
% Подписываем точки
X = 'X%d';
for i=1:length(x1) % цикл от 1 до конца вектора координат
    n = i-1; % Начинаем с X0
    Name = sprintf(X,n); 
    text(x1(i)+0.05,x2(i)+0.05,Name,'Color','red'); % надпись точки
end % конец цикла for