3D plot of nonlinear system of equations in Matlab

Question

I have solved an optimization problem by hand and in Matlab, where I'm using the fsolve function.

There are two unknown variables and one constraint (the two equations in the function is already the first order condition done by hand).

I'm getting a positive solution and a positive value for pi. My problem is that when I plot pi in Matlab for different values of the two variables, I only get negative values of the surface, which does not correspond to the one obtained by fsolve. I.e. I want the maximum value as obtained by fsolve to be on the surface. Any suggestions?

The following is my Matlab code, thank you in advance.

% Calibration
alphaK = 0.5;
alphaL = 0.5;
bbeta = 0.8;
delta = 0.5;
phi = 0.8;
rho = 0.1;
v = 0.001;
w = 1;

guess = [1;1];      %% Start value for solving

sol = fsolve(@function_simple, guess)       %% Solve equations


pi = (alphaK.*sol(1).^rho + alphaL.*sol(2).^rho).^(bbeta./rho) - sol(1) - w.*sol(2)        %% Profit


%% 3D PLOT

% Define a range of values for z(1) and z(2)
z1_range = linspace(0.1, 2, 100);  % Adjust the range as needed
z2_range = linspace(0.1, 2, 100);  % Adjust the range as needed

% Create a meshgrid of z(1) and z(2) values
[z1, z2] = meshgrid(z1_range, z2_range);

% Calculate pi for each combination of z(1) and z(2)
pi = (alphaK .* z1.^rho + alphaL .* z2.^rho).^(bbeta/rho) - z1 - w .* z2;

% Plot the 3D surface plot
figure;
surf(z1, z2, pi);
title('3D Plot of \pi');
xlabel('K');
ylabel('L');
zlabel('\pi');

Function with equations

function F = function_simple(z)

% Calibration
alphaK = 0.5;
alphaL = 0.5;
bbeta = 0.8;
delta = 0.5;
phi = 0.8;
rho = 0.1;
v = 0.001;
w = 1;


% Equations

F(1) = (z(1).*(1+phi.*delta)+z(2).*(phi+w)-v).^((bbeta-rho)./bbeta) - (1-w.*delta)./(alphaK.*bbeta.*z(1).^(rho-1) - delta.*alphaL.*bbeta.*z(2).^(rho-1));

F(2) = alphaK.*z(1).^rho + alphaL.*z(2).^rho - (z(1).*(1+phi.*delta)+z(2).*(phi+w)-v).^(rho./bbeta);

end