'Index in position 1 is invalid. Array indices must be positive integers or logical values.', error in line 14 i.e., t1=D(ax(i),ay(i)).

clc
clear all
syms x y 
f=input('enter the function:')
p=diff(f,x);
q=diff(f,y)
[ax,ay]=solve(p,q)
ax=double(ax); ay=double(ay);
r=diff(p,x); t=diff(q,y); s=diff(p,y);
D=r*t-s^2
figure
fsurf(f)
for i=1:size(ax)
    t1=D(ax(i),ay(i));
    t2=r(ax(i),ay(i));
    t3=f(ax(i),ay(i));
    if double(t1)==0
        sprintf('further investigation required at point:',ax(i),ay(i))
        legstr=[legstr,{'further investigation point'}]
        mkr='ko'
    else
        if double(t1)<0
            sprintf('the saddle point is',ax(i),ay(i))
            sprintf('the value of f at saddle point is',t3)
            legstr=[legstr,{'the saddle point'}]
            mkr='bv'
        else
            if double(t2)<0
                sprintf('the maximum value is at:',ax(i),ay(i))
                sprintf('the maximum value of f :',t3)
                legstr=[legstr,{'point of maximum'}]
                mkr='g+'
            else
                sprintf('the minimum is at',ax(i),ay(i))
                sprintf('the minimum value of f is',t3)
                legstr=[legstr,{'point of minimum'}]
                mkr='r*'
            end
        end
    end
    hold on
    plot3(ax(i),ay(i),t3,mkr)
end
legend(legstr)

[KSSV]

It is asking to input a function. What are the inputs?

[Dyuman Joshi]

Assuming you are giving a function to the input -

size() returns a 2-element vector output, so the loop counter will go from 1 to 1 only.

for i=1:size(ax)

Try -

for i=1:numel(ax)
%or
for i=1:prod(size(ax))

About the error, it seems you are trying to substitute the values in D, r and t to obtain t1, t2 and t3. In that case, I would suggest you check out subs

Or define your functions as functions of x and y, to directly substitute values and get output, for e.g -

syms x y
f(x,y) = sin(x)+cos(y);
f(pi/2,0)

ans = 

2

g(x,y) = sin(x)*cos(y);
g(pi/6,pi/3)

ans = 

h(x,y) = f(x,y)*g(x,y);
h(pi/4,pi/4)

ans = 

Answer 1

[VBBV]

clc
clear all
syms x y 
f=x*y^2;%input('enter the function:')
p=diff(f,x);
q=diff(f,y);
[ax,ay]=solve([p,q],[x y]);
ax=double(ax); ay=double(ay);
r=diff(p,x); 
t=diff(q,y); 
s=diff(p,y);
D=r*t-s^2

D = 

figure
fsurf(f)

The next part of your code seems to be confusing and not sure what you're trying to do