syms x1 x2 

% initiation of  'beta' coeeficient and definition of criterion function

beta = -5 ,  f = x1^2+x2^2+beta*x1*x2 + x1 + 2*x2

%  calculate gradient via directional diferentiation

g = [diff(f,x1);diff(f,x2)]

%  the same using gradient function

g2 = gradient(f,[x1,x2])

% solve zero gradient equations

s = solve([g2(1)==0,g2(2)==0],[x1,x2])
[s.x1, s.x2]

% calculate hessian and its eigenvalues

h = hessian(f,[x1,x2])
lambda = eig(h)

% plot function and critical point

x1 = -5:0.2:5;
x2 = x1';
z = x1.^2+x2.^2+beta*x1.*x2 + x1 + 2*x2;
[px1,px2] = gradient(z);
figure
contour(x1,x2,z,10)
hold on
plot(s.x1,s.x2,'ro')
quiver(x1,x2,px1,px2)
hold off