Maxima Minima

syms x y real
f=input('Enter the function f(x,y):')
dfx=diff(f,x); %first derivative
dfy=diff(f,y);
eqns=[dfx==0,dfy==0];
S=solve(eqns,[x y], 'Real', true);% Gives stationary points
G=S.x;
H=S.y;
dfxx=diff(dfx,x); % f_xx %double derivative
dfyy=diff(dfy,y); % f_yy
dfxy=diff(dfy,x); %f_xy
A=subs(dfxx,S); %for AC-B^2
B=subs(dfxy,S);
C=subs(dfyy,S);
fun=subs(f,S); % Subsitutuion of values
for i=1: length(A)
if A(i)*C(i)-B(i)*B(i)>0 & A(i) < 0
fprintf('The Maximum point is %d %d',[G(i), H(i)]);
fprintf('\n The maximum value is %d',fun(i));
figure;
fsurf(f,[-3 3 -3 3]);
hold on
plot3(G(i),H(i),fun,'*w','markersize',10)
hold off
elseif A(i)*C(i)-B(i)*B(i) > 0 & A(i) > 0
fprintf('Minimum point is %d %d',[G(i), H(i)]);
fprintf('\n The minimum value is %d',fun(i));
figure;
fsurf(f,[-3 3 -3 3]);
hold on
plot3(G(i),H(i),fun,'*w','markersize',10)
hold off
elseif A(i)*C(i)-B(i)*B(i) < 0 fprintf('Saddle point is %d %d', [G(i), H(i)]);
else
fprintf('No conclusion can be drawn for %d %d',[G(i), H(i)]);
end
end