function [x_vals, f_vals, k] = method_SteepDesc_Proj(f, grad_f, xk, sk, limmits, tol, max_iter, mode)
    % f:        Objective function
    % grad_f:   Gradient of the function
    % xk:       Initial point [x0; y0]
    % sk:       Step size (fixed positive scalar)
    % limits:   Bounds of the feasible set for each dimension
    % tol:      Tolerance for stopping criterion
    % max_iter: Maximum number of iterations
    
    % x_vals: Vector with the (x,y) values until minimum
    % f_vals: Vector with f(x,y) values until minimum
    % k:      Number of iterations
    
    if strcmp(mode, 'armijo') == 1
        gamma_f = @(f, grad_f, dk, xk) gamma_armijo(f, grad_f, dk, xk);
    elseif strcmp(mode, 'minimized') == 1
        gamma_f = @(f, grad_f, dk, xk) gamma_minimized(f, grad_f, dk, xk);
    else % mode == 'fixed'
        gamma_f = @(f, grad_f, dk, xk) gamma_fixed(f, grad_f, dk, xk);
    end
    
    % Project the first point if needed
    xk = ProjectionPoint(xk, limmits);

    % Storage for iterations, begin with the first point
    x_vals = xk;
    f_vals = f(xk);
    for k = 1:max_iter
        % Check for convergence
        if norm(grad_f(xk)) < tol
            break;
        end
        dk = - grad_f(xk);          % Steepset descent direction

        % First calculate xk-bar and project it if nessesary
        xkbar = xk + sk * dk;
        xkbar = ProjectionPoint(xkbar, limmits);

        dk = (xkbar - xk);          % Steepest descent projection direction
        
        gk = gamma_f(f, grad_f, dk, xk);   % Calculate gamma
        
        x_next = xk + gk * dk;      % Update step
        f_next = f(x_next);
        
        xk = x_next;                % Update point
        x_vals = [x_vals, x_next];  % Store values
        f_vals = [f_vals, f_next];  % Store function values
    end
end