One problem is using zero for any initial parameter estimate, and especially if the parameter is the only element in the denominator, since that becomes Inf and the solver immediately stops.
Start with different initial estimates instead —
% Define the function f(x1, x2)
f = @(x) x(1)^2 + x(1)*x(2) + 3/2*x(2)^2 - 2*log(x(1)) - log(x(2));
% Define the gradient ∇f(x1, x2)
gradient = @(x) [2*x(1) + x(2) - 2/x(1); x(1) + 3*x(2) - 1/x(2)];
% Define a function that returns a vector for fsolve
stationary_points = fsolve(@(x) gradient(x), rand(2,1));%, options);
% The variable stationary_points now contains the stationary points
disp('Stationary Points:');
disp(stationary_points);
You apparently defined an options structure, however did not include it, so I changed the fsolve call to exclude it.
.
