%% Run an initial plot to get an idea of where to look for minima

clear

xvec = -1:0.01:1;
yvec = xvec;

[XX,YY] = meshgrid(xvec,yvec);
f = 2*XX.^2 - (XX - 1).*(YY - 2).^2 + 2*YY.^4;

mesh(xvec,yvec,f)

%% Newton's Method

clear

%x = [0.1; 0.1];
%x = [1; 0.1];
 x = [1; 0.5];

for n = 1:10000  
    df = [4*x(1)-(x(2)-2)^2 ; -2*(x(1)-1)*(x(2)-2) + 8*(x(2))^3];
    Hf = [4, -2*(x(2)-2) ; -2*(x(2)-2), -2*(x(1)-2) + 24*(x(2))^2];
    
    x = x - Hf\df;
end

fval = 2*(x(1))^2 - (x(1) - 1)*(x(2) - 2)^2 + 2*(x(2))^4;
fprintf(strcat('f= ',num2str(fval)))

%% Calculate critical points analytically

syms xs ys

[xs, ys] = solve([4*xs-(ys-2)^2, -2*(xs-1)*(ys-2) + 8*(ys)^3],[xs,ys]);

xs = double(xs);
ys = double(ys);