syms l
x=sym('x',[1,2])
c=[1 3];
q=[0.3 0.5];
R0=0.98;
Cost=sum(c.*x);
Rel=(1-q(1)^(x(1)+1))*(1-q(2)^(x(2)+1))
L=Cost+l*(Rel-R0)
solve([Rel==R0, gradient(L)==0], [x l])