how to find optimal solutions for the given problem?

I got the answer. but i want to get optimal answer. I couldnt able to get convex graph for this problem.

function newton2D10()
%New trade credit case 1 T>=t
clc
close all
n=0:1:10;
q=0:100:900;
[q,n]=meshgrid(q,n);
TCTC=f(q,n);
mesh(q,n,TCTC)
r0=[0.20;10];
alfa=0.2;
while abs (f(r0(1),r0(2))) > 1e-2
    r0 = r0 - alfa.*(f(r0(1),r0(2)))./fprime(r0(1),r0(2));
end
hold on
r0=[8;0.1]
alfa=0.2;
for index = 1:100
    r0 = r0 - alfa.*inv(fdb1prime(r0(1),r0(2))).*fprime(r0(1),r0(2));
end
r0
f(r0(1),r0(2))
plot(r0(1),r0(2),f(r0(1),r0(2)),'rs','Markersize',20)
    function TC = f(q,n)
        pr=7;
        p=2000;
        P=1;
        Ic=0.15;
        t=0.25;
        Id=0.05;
        %p=1600;
        %pmin=700;
        %pmax=1200;
        %D=1600;
        D=1000;
        z=4;
        m=10;
        hr=5;
        %h1=12;
        %h2=11;
        %h3=13;
        %hb=19.24;
        h1=2;
        h2=1;
        h3=3;
        hb=9;
        cL=4;
        c=2;
        c0=12;
        c1=13;
        c2=12;
        
        ce=15;
        Ct=100;
        k=0.5;
        %k1=0.05;
        k1=70;
        alpha=0.01;
        beta=0.099;
        %u=0.999;
        u=0.999;
        %s=250;
        %A=150;
        s=25;
        A=15;
        g0=5;
        g1=2;
        theta1=1;
        theta2=0.6;
        w0=0.3;
        w1=0.2;
        F=50;
        sigma=5;
        %ts=0.17;
        %pi=80;
        ts=0.25;
        pi=80;
        zeta=0.5;
        %zeta1=300;
        %zeta2=1/300;
        zeta1=300;
        zeta2=1/300;
        sp=15;
        sc=0.5;
        pc=3;
        Ic=0.15;
        t=0.25;
        TC=((((2.*A+2.*F+((D.*P.*t^2).*(Ic-
        Id)))./2)+(n.*Ct)).*(D./n.*q))+hb.*((q./2)+(k1.*sigma.*sqrt(ts+(q./p)))+(((D.
        *pi.*sigma)./(n.*q.*2)).*sqrt(ts+(q./p))).*(sqrt(1+k1^2).*k1))+((zeta+sp).*D)
        + ((n.*q).*((2.*sc+P.*Ic)./2))+(pr.*D)
        +(z.*(m.*D./n.*q))+((hr+cL.*k).*(D.*n.*q./2.*m.*p))+(h1.*(1-(((2.*u-
        1).*D)./(2.*p.*u))).*n.*q)-
        (h1.*n.*q./2)+(s.*D./n.*q)+((D./u).*(c0+(c1.*alpha)+(c0.*alpha)+(c2.*(1-
        beta).*alpha)))+((c.*D./n.*q).*(g0+g1.*n.*q))+(((D.*ce.*n.*theta1)./(n.*q)).*
        (g0+g1.*n.*q))+((w0+w1.*(((1-(((2.*u-1).*D)./(2.*p.*u))).*n.*q)-
        (n.*q./2))).*ce.*theta2)+(D.*zeta1.*p)+(D.*zeta2./p)+((pc.*Ic)./(1+(Ic.*t)))-
        ( D.*P.*t.*Ic);
        
        function TCprime = fprime(q,n)
            dfdq =-((((2*A+2*F+((D*P*t^2)*(Ic-
            Id)))/2)+(n*Ct))*(D/n*q^2))+(hb*(1/2+((k1*sigma)/(2*p*sqrt(ts+(q/p))))))-
            ((D*pi/n*q^2)*((sigma/2)*(sqrt(ts+(q/p))*(sqrt(1+k1^2)-
            k1))))+((D*pi/n*q)*((sigma/2*p*sqrt(ts+(q/p))*(sqrt(1+k1^2)-
            k1))))+(n*((2*sc+P*Ic)/2))-(z*m*D/n*q^2)+((hr+cL*k)*(D*n/2*m*p))+(h1*(1-
            
            (((2*u-1)*D)/(2*p*u)))*n)-(h1*n/2)-(s*D/n*q^2)-(c*D*g0/n*q^2)-
            (D*ce*n*theta1*g0/n*q^2)+(w1*(1-(((2*u-1)*D)/(2*p*u))*n-(n/2))*ce*theta2);
            dfdn =-((((2*A+2*F+((D*P*t^2)*(Ic-Id)))/2))*(D/q*n^2))-
            ((D*pi/q*n^2)*((sigma/2)*(sqrt(ts+(q/p))*(sqrt(1+k1^2)-
            k1)))+(q*((2*sc+P*Ic)/2))-(z*m*D/q*n^2)+((hr+cL*k)*(D*q/2*m*p))+(h1*(1-
            (((2*u-1)*D)/(2*p*u)))*q)-(h1*q/2)-(s*D/q*n^2)-
            (c*D*g0/q*n^2)+D*ce*theta1*g1)+(w1*(1-(((2*u-1)*D)/(2*p*u))*q-
            (q/2))*ce*theta2);
            
            TCprime =[dfdq;dfdn]
            function TCdb1prime = fdb1prime(q,n)
                df2dq2=((((2*A+2*F+((D*P*t^2)*(Ic-
                Id)))/2)+n*Ct)*(2*D/n*q^3))+((hb*k1*sigma/2*p)*((ts+(q/p)^(-3/2))))
                +((2*D*pi/n*q^3)*((sigma/2)*(sqrt(ts+(q/p))*(sqrt(1+k1^2)-
                k1))))+((D*pi/n*q^2)*(sigma/(2*p*sqrt(ts+(q/p))))*(k1-
                sqrt(1+k1^2)))+((D*pi/n*q^2)*(sigma/2*p*sqrt(ts+(q/p)))*(k1-
                sqrt(1+k1^2)))+((D*pi/n*q)*(sigma/4*p^2)*((ts+(q/p))^(-3/2))*(k1-
                sqrt(1+k1^2)))+(2*z*m*D/n*q^3)+(2*s*D/n*q^3)+(2*c*D*g0/n*q^3)+(2*D*ce*theta1*
                g0/q^3);
                df2dqn=((((2*A+2*F+((D*P*t^2)*(Ic-
                Id)))/2))*(2*D/q^2*n^2))+((2*sc+P*Ic)/2)-
                ((D*pi/n^2*q^2)*((sigma/2)*(sqrt(ts+(q/p))*(sqrt(1+k1^2)-
                k1))))+((D*pi/n^2*q)*(( sigma/2*p*sqrt(ts+(q/p))*(sqrt(1+k1^2)-k1))))+sc-
                (z*m*D/n^2* q^2)+((hr+cL*k)*(D/2*m*p))+(h1*(1-(((2*u-1)*D)/(2*p*u))))-(h1/2)-
                (s*D/n^2*q^2)-(c*D*g0/n^2*q^2)+w1*(1-(((2*u-1)*D)/(2*p*u))-(1/2))*ce*theta2;
                df2dnq =((((2*A+2*F+((D*P*t^2)*(Ic-
                Id)))/2))*(2*D/q^2*n^2))+((2*sc+P*Ic)/2)-
                ((D*pi/n^2*q^2)*((sigma/2)*(sqrt(ts+(q/p))*(sqrt(1+k1^2)-
                k1))))+((D*pi/n^2*q)*(( sigma/2*p*sqrt(ts+(q/p))*(sqrt(1+k1^2)-k1))))+sc-
                (z*m*D/n^2* q^2)+((hr+cL*k)*(D/2*m*p))+(h1*(1-(((2*u-1)*D)/(2*p*u))))-(h1/2)-
                (s*D/n^2*q^2)-(c*D*g0/n^2*q^2)+w1*(1-(((2*u-1)*D)/(2*p*u))-(1/2))*ce*theta2;
                df2dn2=((((2*A+2*F+((D*P*t^2)*(Ic-Id)))/2))*(2*D/q*n^3))-
                ((D*pi/n^2*q^2)*((sigma/2)*(sqrt(ts+(q/p))*(sqrt(1+k1^2)-
                k1))))+((D*pi/n^2*q)*(( sigma/2*p*sqrt(ts+(q/p))*(sqrt(1+k1^2)-k1))))+sc-
                (z*m*D/n^2* q^2)+((hr+cL*k)*(D/2*m*p))+(h1*(1-(((2*u-1)*D)/(2*p*u))))-(h1/2)-
                (s*D/n^2*q^2)-(c*D*g0/n^2*q^2)+w1*(1-(((2*u-1)*D)/(2*p*u))-(1/2))*ce*theta2;
                
            end
        end
    end
end