function [C]=C_steady_state(alpha, delta, p_over_p_N, iota, kappa, gamma, habit, beta, mc_N, p_over_p_T)

C = fzero(@(C) (1-alpha)*((delta*C*p_over_p_N^iota)^kappa)*(alpha*(((1-delta)*C*(p_over_p_T^iota))/(delta*C*(p_over_p_N^iota)))^gamma...
    + 1 - alpha)^((1+kappa-gamma)/gamma)- ((1-habit*beta)/((1-habit)*C))*mc_N/(p_over_p_N) , [0.1, 10]); 