function [xopt,fopt,niter,gnorm,dx] = grad_descent
% min f(x)=x1^2+x1x2+3x2^2
%%%%%%%%%Initial value%%%%%%
x0=[3,3]';
% termination tolerance
tol = 1e-8;
% maximum number of allowed iterations
maxiter = 10000;
% minimum allowed perturbation
dxmin = 1e-8;

% initialize gradient norm, optimization vector, iteration counter, perturbation
gnorm = inf; 
x = x0; 
niter = 0; 
dx = inf;


% define the objective function:
f = @(x1,x2) x1.^2 + x1.*x2 + 3*x2.^2;


% plot objective function contours for visualization:
[X,Y]=meshgrid(-5:0.1:5,-5:0.1:5);
figure(1)
contour(X,Y,f(X,Y),50,LineWidth=1)
hold on


% redefine objective function syntax for use with optimization:
f2 = @(x) f(x(1),x(2));


% gradient descent algorithm:
while (gnorm>=tol&&niter <= maxiter&&dx >= dxmin)
    % calculate gradient:
    g = grad(x);
    gnorm = norm(g);

    %%%%%%%Line search%%%%%%%%
    t=1;
    xnew=x-t*g;
     while f2(xnew)>f2(x)-0.9*t*gnorm^2
         t=0.2*t;
         xnew=x-t*g;
     end

%     %%Constant stepsize:%%%%%
%       alpha = 0.2;
%       xnew = x - alpha*g;

    % plot current point
     plot([x(1) xnew(1)],[x(2) xnew(2)],'ko-')
     
    % update termination metrics
    niter = niter + 1;
    dx = norm(xnew-x);
    x = xnew;
    
end
hold off
xopt = x;
fopt = f2(xopt);
norm(g);
niter = niter - 1;
% define the gradient of the objective
function g = grad(x)
 g = [2*x(1) + x(2);
     x(1) + 6*x(2)];
end
end