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

psi1  = upar(1);
psi2  = upar(2);
psi3  = upar(3);
sigma = upar(4);

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

