function [F] = myfunc (x) global beta gamma P y Uau Upr w=x(1:2,1); c=x(3:4,1); uc=(c.^(1-gamma))/(1-gamma); F(1:2,1) = uc + beta*P*w - w; F(3,1) = [1 1]*c - 1; F(4,1) = w(2) - max(Uau(2),Upr);