%% Singularity of Hessian
% The task is to minimize the function
%  f(x_1,x_2) = x_1^2 + x_2^4
%% Computing the stationary point

syms x1 x2

f = x1^2 + x2^4

g = gradient(f)

sol = solve(g,[x1,x2])

xst = [sol.x1; sol.x2]

%
%% Computing the Hessian

H = hessian(f)

%
% and evaluating it at the stationary point

Hst = subs(H,[x1;x2],xst)

%
%% and higher-order derivatives of the original function

fxxx(:,:,1) = diff(H,x1)
fxxx(:,:,2) = diff(H,x2)

fxxxx(:,:,:,1) = diff(fxxx,x1)
fxxxx(:,:,:,2) = diff(fxxx,x2)

%
% Obviously the fourth partial derivative \frac{\partial^4 f(x_1,x_2)}{\partial x_2^4} is positive, hence the stationary point is a minimum.
%% Plotting

x1 = -10:.01:10;
x2 = x1;

[X1,X2] = meshgrid(x1,x2);

F = X1.^2 + X2.^4;

[C,h] = contour(X1,X2,F);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
xlabel('x_1'), ylabel('x_2')
grid
hold on
plot(xst(1),xst(2),'ro')
hold off