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

psi1  = upar(1);
psi2  = upar(2);
psi3  = upar(3);
sigma = upar(4);
u     = (cm.^psi1.*cn^psi2.*lei^psi3).^(1-sigma)/(1-sigma);
u1    = psi1*(1-sigma)*u./cm;
u2    = psi2*(1-sigma)*u./cn;
u3    = psi3*(1-sigma)*u./lei;

