syms x y;
f = (1 - x)^2 + 100 * (y - x^2)^2;

[x, xs, gs] = gd(f, [0;0], 500, 1e-6, true);
display(x);
display(xs);
display(gs);
plot(xs(1, :), xs(2, :), '-o');

function [xk, xs, gs] = gd(f, x0, max_itr, tol, bt)
    syms x y;
    xk = x0;
    gf = gradient(f);
    hf = hessian(f, [x, y]);
    tnorm = norm(at(gf, xk), 2);
    xs = zeros(size(x0, 1), max_itr + 1);
    xs(:, 1) = x0;
    gs = zeros(1, max_itr + 1);
    gs(1) = tnorm;
    k = 1;
    
    while tnorm > tol && k <= max_itr 
        p = -at(gf, xk);
        H = at(hf, xk);
        if bt
            alpha = backtracking(f, xk, 1, p, 0.01, 1e-4)
        else
            alpha = dot(-H * xk, p) / dot(H * p, p)
        end
        xk = xk + alpha * p
        k = k + 1;
        xs(:, k) = xk;
        gs(k) = tnorm;
    end
    xs = xs(:, 1:k);
    gs = gs(:, 1:k);
end