clear all
close all


set(0, 'DefaultLineLineWidth', 1, ...
    'DefaultLineMarkerSize', 6,...
    'DefaultAxesFontSize', 18, ...
    'DefaultAxesTickLabelInterpreter', 'latex',...
    'DefaultFigurePaperOrientation', 'portrait',...
    'DefaultFigurePaperUnits', 'inches',...
    'DefaultFigurePaperSize', [6, 4.5],...
    'DefaultFigurePaperPosition', [0,0,6,4.5],...
    'DefaultLegendInterpreter', 'latex',...
    'DefaultTextInterpreter', 'latex',...
    'DefaultFigureRenderer', 'painters')

% Set up the grid
[x, y] = meshgrid(-3:.05:3, -3:.05:3);  % Tighter domain to match radius ~2.24
x = x';
y = y';

% Compute the objective function
z = (x.^2) .* exp(y) + 1.0 .* x .* (y.^2);

% Define custom contour levels
custom_levels = [-32 -16 -8 -4 -2 -1 1 4 8 16 25 48];

% Plot contours
close all;
figure;
set(gcf, 'Position', [10, 10, 800, 800]);
 
subplot(1,2,1)
colormap turbo
contour(x, y, z, custom_levels, 'ShowText', 'on', 'LineWidth', 3);
hold on;

% Plot constraint: x^2 + y^2 = 5
theta = linspace(0, 2*pi, 300);
r = sqrt(5);
xc = r * cos(theta);
yc = r * sin(theta);
zc = (xc.^2).* exp(yc) + 1.0*xc.*(yc.^2);
plot(xc, yc, 'k', 'LineWidth', 2, 'LineStyle','--');

xline(0, 'k');
yline(0, 'k');
xlabel('$x$')
ylabel('$y$')
grid
axis equal;
title('Top View', 'Interpreter','latex', 'FontSize', 20);

% Plot optimal points as stars
[~, maxIdx] = max(zc);
x_opt = xc(maxIdx);
y_opt = yc(maxIdx);
z_opt = zc(maxIdx);
plot(x_opt, y_opt, 'kp', 'MarkerSize', 12, 'MarkerFaceColor', 'k');  % Star at (-x, y)

%
subplot(1,2,2) 
% Define grid
[x, y] = meshgrid(-2.2:0.25:2.2, -2.2:0.25:2.2);
F = (x.^2).* exp(y) + 1.0*x.*(y.^2);

% Plot the surface
meshc(x, y, F);
colormap turbo
hold on;
set(gcf, 'Position', [46 107 767 682]);

% Now plot the constraint curve g(x, y) = 5
theta = linspace(0, 2*pi, 400);
r = sqrt(5);
xc = r * cos(theta);
yc = r * sin(theta);
zc = (xc.^2).* exp(yc) + 1.0*xc.*(yc.^2);

plot3(xc, yc, zc+.1, 'k', 'LineWidth', 2);  % Thick black constraint curve on surface
plot3(xc, yc, -30*ones(1,400), 'k--', 'LineWidth', 1.5);% Labels
xlabel('$x$'); ylabel('$y$'); zlabel('$F(x,y)$');
title('$F(x,y) = x^2 e^y$ with constraint curve $x^2 + y^2 = 5$');
view(-10, 10);  % nice 3D angle
 title('Side View', 'Interpreter','latex', 'FontSize', 20);


% Find maximum on the constraint
[~, maxIdx] = max(zc);
x_opt = xc(maxIdx);
y_opt = yc(maxIdx);
z_opt = zc(maxIdx);

% Plot optimal point as red star
plot3(x_opt,  y_opt, z_opt, 'kp', 'MarkerSize', 12, 'LineWidth', 2);
plot3(x_opt,  y_opt, -30, 'kp',   'MarkerSize', 12, 'LineWidth', 2);

    x0=10;
    y0=10;
    width=1000;
    height=500;
    set(gcf,'position',[x0,y0,width,height])
    tightfig
saveas(gcf,'lect2_lagrange2.pdf')

%% THE END