function [u,u1,u2,u3] = util2(cm,cn,lei,upar);

psi1  = upar(1);
psi2  = upar(2)+1e-10;
psi3  = upar(3);
sigma = upar(4);

alp   = 1;
C     = (psi1*cm.^psi2+(1-psi1)*cn.^(psi2*alp)).^(1/psi2);
u     = (C.^(1-psi3).*lei.^psi3).^(1-sigma)/(1-sigma);
u1    = (1-psi3)*(1-sigma)*psi1*u.*cm.^(psi2-1)./ ...
          (psi1*cm.^psi2+(1-psi1)*cn.^(psi2*alp));
u2    = alp*(1-psi3)*(1-sigma)*(1-psi1)*u.*cn.^(psi2*alp-1)./ ...
          (psi1*cm.^psi2+(1-psi1)*cn.^(psi2*alp));
u3    = psi3*(1-sigma)*u./lei;

