function x = fxsteep(f,e,a,b)
x1 = a;
x2 = b;
Q = fxhesson(f,x1,x2);
x0 = [x1,x2]';
temp = [x0];
fx1 = diff(f,'x1');
fx2 = diff(f,'x2');
g = [fx1,fx2]';
g1 = subs(g);
d = - g1;
while (abs(norm(g1))> e)
lamda = (-d)'*d/((-d)'*Q*d);
x0 = x0 - lamda *g1;
temp = [temp;x0];
a = [1,0] * x0;
b = [0,1] * x0;
x1 = a;
x2 = b;
g1 = subs(g);
d = - g1;
end
x = temp;
end

function Q = fxhesson(f,a,b)
x1 = a;
x2 = b;
fx1 = diff(f,1,'x1');
fx2 = diff(f,1,'x2');
fx1x1 = diff(fx1,1,'x1');
fx1x2 = diff(fx1,1,'x2');
fx2x2 = diff(fx2,1,'x2');
fx2x1 = diff(fx2,1,'x1');
fx1x1 = subs(fx1x1);
fx1x2 = subs(fx1x2);
fx2x1 = subs(fx2x1);
fx2x2 = subs(fx2x2);
Q = [fx1x1,fx1x2;fx2x1,fx2x2];

end

syms x1 x2;
X=[x1,x2];
fx=X(1)^2+4*X(2)^2;
z=fxsteep(fx,0.0001,1,1)

ans

1
1
48/65
-3/65
36/325
36/325
1728/21125
-108/21125
1296/105625
1296/105625
62208/6865625
-3888/6865625
46656/34328125
46656/34328125
2239488/2231328125
-139968/2231328125
1679616/11156640625
1679616/11156640625
80621568/725181640625
-5038848/725181640625
60466176/3625908203125
60466176/3625908203125
2902376448/235684033203125
-181398528/235684033203125