I want to find two unknown variables using newton function method

I couldn't find q and n value from the first derivative of this coding. Kindly do the needful as soon as possible.

function Hemamam()
clc
close all
q=100:10:550
n=5:10:175
[qq,nn]=meshgrid(q,n);
TCTC=f(qq,nn);
mesh(qq,nn,TCTC)
r0=[125;15];
alfa=0.5;
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
    function TC = f(q,n)
        pr=7;
        p=2500;
        D=1600;
        z=4;
        m=10;
        hr=5;
        h1=12;
        h2=11;
        h3=13;
        cL=4;
        c=20;
        c0=2;
        c1=3;
        c2=2;
        ce=15;
        k=0.1;
        alpha=0.1;
        beta=0.99;
        u=0.999;
        s=250;
        A=110;
        g0=15;
        g1=20;
        theta1=1;
        theta2=0.6;
        w0=0.3;
        w1=0.2;
        TC=(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.*A./n.*q)+(D.*h3./q)+(D.*h2./2.*n)+((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)
        function TCprime = fprime(q,n)
        end
    end
P1=diff(TC,q)
P2=diff(TC,n)
dfdq=pretty(P1)
dfdn=pretty(P2)
TCprime =[dfdq;dfdn];
syms q
q=solve(P1,q);
q
syms n
n=solve(P2,n);
n
end