function F = steadystate(x, csistar, etastar, Ystar, csi, eta, omegabar, chi, psi, beta, s, alpha, tauf, sigman, sigmau, etal, tauw, tauc, epsilon, A, gi, rho, delta, ti);

F = [
% 1
-((1 - csistar)*(csistar*x(1)^(1 - etastar) + (1 - csistar))^(1/(1 - etastar))*Ystar) + (1 - csi)*x(1)^(1 - eta)/(csi + (1 - csi)*x(1)^(1 - eta))*x(2);
% 2
-(x(3)*omegabar/chi*x(4)^(1 - psi)) + beta/(1 - beta*(1 - s))*x(2)/x(5)*((1 - alpha) - (1 + tauf)*x(6)*x(7)*x(5)/x(2));
% 3
-(sigman*x(7)^(1 + etal)) + (1 - tauw)/((1 + tauc)*(1 + tauf))*(1 - alpha)/x(5)*(1/(x(3)*x(8)/x(2)));
% 4
-(sigmau*x(9)^(1 + etal)) + (1 - tauw)/((1 + tauc)*(1 + tauf))*(1 - epsilon)/epsilon*x(3)*omegabar*x(9)*x(4)/x(2)*(1/(x(3)*x(8)/x(2)));
% 5
-x(14) + (alpha*beta/(1 - beta*(1 - delta)))^(1/(1 - alpha))*A^(1/(1 - alpha))*x(3)^(-1/(1 - alpha));
% 6
-x(3) + (csi + (1 - csi)*x(1)^(1 - eta))^(1/(1 - eta));
% 7
-gi + tauc*(x(3)*x(8)/x(2)) + (tauf + tauw)*x(6)*x(7)*x(5) - (1 - tauw)*rho*x(6)*x(7)*(1 - x(5))/x(2) - ti;
% 8
-x(2) + A*(x(14))^alpha*x(5)*x(7);
% 9
-(x(6)*x(7)*x(5)/x(2)) + (etal*(1 - epsilon)*x(3)*omegabar*x(9)*x(4)*x(5)/x(2) + (1 + etal*(1 - epsilon))*(1 - alpha))/((1 + tauf)*(1 + etal)*(1 - rho*x(9)));
% 10
-s*x(5) + x(9)*x(10)*(1 - x(5));
% 11
-x(11) + chi*(x(10)/x(11))^(psi - 1);
% 12
-1 + x(3)*x(8)/x(2) + x(3)*x(12)/x(2) + x(3)*omegabar*x(13)/x(2) + gi;
% 13
-(x(3)*x(12)/x(2)) + alpha*beta*delta/(1 - beta*(1 - delta));
% 14
-x(13) + x(4)*x(9)*(1 - x(5));
% 15
-x(14) + x(15)/(x(5)*x(7))];



% 1
%(1 - csistar)*(csistar*Pf^(1 - etastar) + (1 - csistar))^(1/(1 - etastar))*Ystar = (1 - csi)*Pf^(1 - eta)/(csi + (1 - csi)*Pf^(1 - eta))*Y;
% 2
%P*omegabar/chi*theta^(1 - psi) = beta/(1 - beta*(1 - s))*Y/N*((1 - alpha) - (1 + tauf)*w*h*N/Y);
% 3
%sigman*h^(1 + etal) = (1 - tauw)/((1 + tauc)*(1 + tauf))*(1 - alpha)/N*(1/(P*C/Y));
% 4
%sigmau*e^(1 + etal) = (1 - tauw)/((1 + tauc)*(1 + tauf))*(1 - epsilon)/epsilon*P*omegabar*e*theta/Y*(1/(P*C/Y));
% 5
%k = (alpha*beta/(1 - beta*(1 - delta)))^(1/(1 - alpha))*A^(1/(1 - alpha))*P^(-1/(1 - alpha));
% 6
%P = (csi + (1 - csi)*Pf^(1 - eta))^(1/(1 - eta));
% 7
%gi = tauc*(P*C/Y) + (tauf + tauw)*w*h*N - (1 - tauw)*rho*w*h*(1 - N)/Y - ti;
% 8
%Y = A*(k)^alpha*N*h;
% 9
%w*h*N/Y = (etal*(1 - epsilon)*P*omegabar*e*theta*N/Y + (1 + etal*(1 - epsilon))*(1 - alpha))/((1 + tauf)*(1 + etal)*(1 - rho*e));
% 10
%s*N = e*p*(1 - N);
% 11
%p = theta*q;
% 12
%q = chi*theta^(psi - 1);
% 13
%1 = P*C/Y + P*I/Y + P*omegabar*V/Y + gi;
% 14
%P*I/Y = alpha*beta*delta/(1 - beta*(1 - delta));
% 15
%V = theta*e*(1 - N);