function r = foc(n1,param);

l     = length(n1);
k1    = param(:,1);
k2    = param(:,2);
z     = param(:,3);
tauc  = param(:,4);
tau1  = param(:,5);
tau2  = param(:,6);
taun  = param(:,7);
g     = param(:,8);
thet1 = param(:,9);
thet2 = param(:,10);
x     = param(:,11);
a1    = param(1,12);
a2    = param(1,13);
rho   = param(1,14);
psi   = param(1,15);

n2    = ( a2/a1* (1-thet2)./(1-thet1).*((k2./z).^(thet2*rho)./ ...
          ((k1./z).^(thet1*rho))).*n1.^(1-rho*(1-thet1)) ).^ ...
          (ones(l,1)./(1-rho*(1-thet2)));

y1    = k1.^thet1.*(z.*n1).^(1-thet1);
y2    = k2.^thet2.*(z.*n2).^(1-thet2);
y     = (a1*y1.^rho+a2*y2.^rho).^(1/rho);
c     = y-x-g;
p1    = a1*(y./y1).^(1-rho);

r     = (1-n1-n2).*(1-taun).*p1.*(1-thet1).*(k1./n1).^ ...
             thet1/psi.* z.^(1-thet1)./(1+tauc)-c;

