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);