function  u = cumcu(x,p);

alpha = p(1);
beta  = p(2);
delta = p(3);
gamma = p(4);
lam   = p(5);
omega = p(6);
pa    = p(7);
pb    = p(8);
psi   = p(9);
sig   = p(10);
theta = p(11);
cg    = p(13);
army  = p(14);
k     = x(1);
nl    = x(2);
lz    = x(3);
c1    = x(5);
i     = x(6);
n     = x(7);
h     = x(8);

d  = (n-pa)*pb;
pn = (d*d*d+abs(d*d*d));
c0 = (exp(lz)*lam*k^theta*n^(1-theta)*h-cg-i-c1*n-alpha*(n-nl)^2)/(1-n);
u  = (n*(c1^gamma*(1-h)^(1-gamma))^(1-omega)/(1-omega)+(1-n)*c0^ ...
     (gamma*(1-omega))/(1-omega)-pn)*(1-army);

