how to find multi variables value for nonlinear equation

hi,

I want to find N,Q,and K1 value using n,q,k1 equations using iteration process.  After finding, N,Q and K1 substitute these values in TC, I want to find TC value also.  I dont know where i left mistake in this program.  anyone find the answer and let me know.  
we may consider n, q and k1 any value. after iteration only, we can find N, Q and k1.  n=partial derivative with respect to TC n
clc
clear all
pr=7;
p=700;
%p=2500;
%pmin=700;
%pmax=1200;
D=1600;
z=4;
m=10;
hr=5;
h1=12;
h2=11;
h3=13;
hb=19.24;
cL=4;
c=20;
c0=2;
c1=3;
c2=2;
ce=15;
Ct=100;
k=0.01;
alpha=0.1;
beta=0.99;
u=0.999;
s=250;
A=150;
g0=15;
g1=20;
theta1=1;
theta2=0.6;
w0=0.3;
w1=0.2;
F=50;
sigma=5;
ts=0.17;
pi=80;
zeta=0.5;
sp=15;
sc=0.5;
TC1=inf
j=1;
l=1;
 for i=1:1:20
q=1:10:191
n=1:1:20
k1=1:1:20
k1=sigma.*(sqrt(ts+(q./p)).*(hb+((D.*pi./n.*q.*2).*((k1./sqrt(1+k1.^2))-1))));
q=-((A+F+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.*sc)-(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;
n=-((A+F).*(D./q.*n.^2))-((D.*pi./q.*n.^2).*((sigma./2).*(sqrt(ts+(q./p)).*(sqrt(1+k1.^2)-k1)))+(q.*sc)-(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);
%n=sqrt((D.*zr.*m+s+c.*g0+(A.*((q.*h2)/2))))/(q.*(((hr+cL.*k).*((D.*q)/(m.*p)))+((h1.*q).*((1/2)-(((2.*u-1).*D)/(2.*p.*u))-D))+(D.*ce.*g1.*theta1)+(ce.*((w1.*(1-(((2.*u-1).*D)./(2.*D.*u))).*q)-(q./2)).*theta2)));
%q=sqrt((2.*D.*(zr.*m+s+g0.*(c+n.*ce.*theta1)+A+n.*h3))./(n.^2.*(((hr+cL.*k).*(D./m.*p))+(h1.*(1-(((2.*u-1).*D)./(2.*p.*u)))+(ce.*(w1.*(2-(((2.*u-1).*D)./(p.*u))-1).*theta2))))));
K1(i+1)=k1;
N(i+1)=n;
Q(i+1)=q;
j=j+1;
 end
K1
N
Q
z1=1;
for z=2:1:20
      if ((N(z)~=N(z-1)) | Q(z)~=Q(z-1) | K1(z)~=K1(z-1))
                   z1=z1+1;
               else
                   break
 end
end
z1;
j=1;
for i=1:20
%TC1=(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(i))+((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;
TC1=(A+F+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.*sc)+(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)
TC(j)=TC1
j=j+1;
end
round(TC)
for m=2:20
if(TC1 >=TC)
   TC1=TC;
else
    ansn = m+1;
   break
end
end
round(TC)