To

clc;
clear all;

x1=10;
x2=20;

T=10^(-6);

eps=10^-8;

xk=[x1; x2];

rk=[ (f(x1+T,x2)-f(x1,x2))/T ;(f(x1,x2+T)-f(x1,x2))/T ]

x1
x2


while norm(rk)>eps

hk1= (f(x1+2*T,x2)-2*f(x1+T,x2)+f(x1,x2))/T^2;
hk2= (f(x1+T,x2+T)-f(x1+T,x2)-f(x1,x2+T)+f(x1,x2))/T^2;
hk3= (f(x1+T,x2+T)-f(x1+T,x2)-f(x1,x2+T)+f(x1,x2))/T^2;
hk4= (f(x1,x2+2*T)-2*f(x1,x2+T)+f(x1,x2))/T^2;

H=[hk1 hk2; hk3 hk4]
theta=(rk'*rk)/(rk'*H*rk);

xk=xk-theta*rk;
x1=xk(1);
x2=xk(2);
rk=[ (f(x1+T,x2)-f(x1,x2))/T ;(f(x1,x2+T)-f(x1,x2))/T ]

end

%theta=(rk'*rk)/(rk'*(H*rk));

%q=xk+theta*hk
%q
xk