function k=steady_state(initial)

c=initial(1);
ef=initial(2);
em=initial(3);
initial


% parameters

rho = 0.97;
theta = 0.4;
kappa = 0.33;
gammaf = 0.6;
gammam = 0.8;
alpha = 0.3;
beta = 0.3;
nu = 0.15;
eta = 0.3;
phi = 1;

r= 1/rho-1;
k = (r/kappa)^(1/(kappa-1));
wf = eta*(1-kappa)*k^kappa;
wm = (1-eta)*(1-kappa)*k^kappa;
n = 2*(ef/em)^(1/(gammaf-gammam));


k(1) = c-k*(1+r)*(1 + (ef/em)^(1/(gammaf-gammam)));
k(2) = (gammaf*phi/2)*(ef/em)^((gammaf-1)/(gammaf-gammam))*((1-beta)*log(ef) + beta*log(em)) + (gammam*phi/2)*(ef/em)^((gammam-1)/(gammaf-gammam))*(alpha*log(ef) + (1-alpha)*log(em)) - (1/c)*(nu*wf*ef^(1-beta)*em^beta + 0.5*ef + 0.5*em);
k(3) = c + (ef/em)^(1/(gammaf-gammam))*ef + (ef/em)^(1/(gammaf-gammam))*em + c/(1+r) - wf*ef^(1-beta)*em^beta*(1-nu*2*(ef/em)^(1/(gammaf-gammam))) - wm*ef^alpha*em^(1-alpha);

end