why solve function gives three results instead of one

I ran the code below and try to solve a Lagrangian multiplier problem. I hope to get one result from this problem, but I unexpectedly get three results instead of one. Please help me check what potential problem could be.

* Setup parameters;
lambda=0.5;
Rbar=2;
Rb=1.4;
ro=2;

syms y gamma
e1 = 2*(y/lambda)^(-ro)-((1-y)/(1-lambda)*Rb)^(-ro)*Rb ...
  -((1-y)/(1-lambda)*(2*Rbar-Rb))^(-ro)*(2*Rbar-Rb)...
  +gamma*((Rb*lambda+1-lambda))==0;
e2 = gamma*((Rb*lambda+1-lambda)*y-Rb*lambda)==0;

sol = solve([e1, e2], [y, gamma]);
ySol = sol.y;
gammaSol = sol.gamma;

Res=[double(ySol) double(gammaSol)];

Answer 1

[Walter Roberson]

The case where gamma = 0 has two solutions, and there is additional solution for non-zero gamma.

[y = 7/12, gamma = 60/637]
[y = 91/41-(5/41)*sqrt(182), gamma = 0]
[y = 91/41+(5/41)*sqrt(182), gamma = 0]

Your second equation is of the form gamma * something = 0 . That is going to have solutions when gamma = 0, and it is going to have solutions when the something = 0, so you should be expecting multiple solutions.

[Xin CUI]

Thanks Walter. That makes sense. But I would like to get one optimal solution. Could you suggest how to sort out which one is the true solution?

Actually, I have a utility function, and would like to maximize it. The utility function is

 u = 2*(y/lambda)^(1-ro)/(1-ro)+(1-lambda)*((1-y)/(1-lambda)*Rb)^(1-ro)/(1-ro)+(1-lambda)*((1-y)/(1-lambda)*(2*Rbar-Rb))^(1-ro)/(1-ro)

And subject to

(Rb*lambda+1-lambda)*y-Rb*lambda<0

I would like to choose y to maximize u,

I can find correct results by using fmincon. Now I would like to use Lagrangian multiplier and write out first order order to solve the problem again. But now I got three solutions.

[Walter Roberson]

Your (Rb*lambda+1-lambda)*y-Rb*lambda<0 has only y as the free variable, so it can be re-arranged as a simple upper bound on y of 7/12 .

When you substitute the constants into u and simplify, you end up with (1/91)*(91-66*y)/(y*(-1+y)) . This has singularities at 0 and 1 and is negative between those two and goes to +infinity at the singularities. The upper bound on y falls between the two singularities. Therefore the maximum u is y as close as you can get to 0 from below.

The same situation holds for a broad range of values, changing only if lambda >= 1 or Rb >= 2*Rbar or ro goes negative. As long as those broad constraints are met, u is maximized where y is the negative number as close as possible to 0.