for n = .0001:.0001:10
M1_o = n;
for k = .0001:.0001:10
M2_o = k;
for r = .0001:.0001:10
N2_o = r;
for iter = 1: itermax
if (iter <= itermax) && (dist > crit1)
M(1) = M1_o;
M(2) = M2_o;
NS(2) = N2_o;
Y(1) = cu(1)^((1-alpha)*mu) * sKM(1)^((1-alpha)*mu) * NS(1)^(alpha*mu) * M(1)^(1-alpha*mu);
Y(2) = cu(2)^((1-alpha)*mu) * sKM(2)^((1-alpha)*mu) * NS(2)^(alpha*mu) * M(2)^(1-alpha*mu);
P(1) = M(1)/((1-mu)*Y(1));
P(2) = (1 - (omega(1)^(1/psi))*(P(1)^((psi-1)/psi)))^(psi/(psi-1))/((1+taxd.d(1))*(1-omega(1))^(1/(psi-1)));
C(1) = sCM(1)*M(1);
X(1) = sXM(1)*M(1);
K(1) = sKM(1)*M(1);
YC(1) = sYM(1)*M(1);
B(1) = sBM(1)*M(1);
y11 = omega(1)^(1/psi) * P(1)^(-1/psi) * YC(1);
y12 = (1-omega(1))^(1/psi) * (P(2)*(1+taxd.d(1)))^(-1/psi) * YC(1);
y21 = Y(1) - y11;
if y21 < 0
warning('y21 is negative')
end
y22 = Y(2) - y12;
if y22 < 0
warning('y22 is negative')
end
YC(2) = ((1-omega(2))*y21^(1-psi) + omega(2)*y22^(1-psi))^(1/(1-psi));
PC(2) = ((1-omega(2))^(1/psi)*(P(1)*(1+taxd.d(2)))^((psi-1)/psi) + omega(2)^(1/psi) * P(2)^((psi-1)/psi) )^(psi/(psi-1));
C(2) = (M(2)/a(2)) * ((1-NS(2))/NS(2)) * (mu*alpha/(1-mu)) * ((1-taxd.n(2))/(1+taxd.c(2)));
sCM(2) = C(2)/M(2);
sBM(2) = (1/(1-mu) - sCM(2) - sXM(2) - 1 - sGM(1) )*PC(2)/(gammax*q -1);
sB(2) = (1-mu)*sBM(2);
sB(1) = -sB(2);
M1_n = (C(1)+X(1)+sG(1,1)+ (gammax*q-1)*sB(1))*((mu-1)/mu);
M2_n = (YC(2) + (gammax*q-1)*sB(2)*(1/PC(2)))*(1-mu);
N2_n = (cu(2)^((1-alpha)*mu) * (sKM(2))^((1-alpha)*mu)* M(2)^(1-alpha*mu) * Y(2)^(-1) )^(-1/(alpha*mu));
dist = max(max(abs(M1_o - M1_n), abs(M2_o - M2_n)), abs(N2_o-N2_n));
M1_o = 0.8 .* M(1) + 0.2 .* M1_n;
M2_o = 0.8 .* M(2) + 0.2 .* M2_n;
N2_o = 0.8 .* N2_o + 0.2 .* N2_n;
end
end
if y21 > 0 && y22 > 0
disp(n)
disp(k)
disp(r)
break
end
end
end
end