@Carl Love @Kitonum

Thank you so much. I used the upward green arrow to attach the file. It seems I made a mistake somewhere. I send the PDF and .mw files both. Thank you so much.
 

Download Code_JN_major.mw

@Kitonum Thanks a million. Your example is exactly what I need. But I received errors when running P3, P4, P5, and P6. I have attached my code. Would you please kindly have a look? 

@acer 

Thank you. Sorry for this (I did not know I shouldn't edit the question significantly).

Also, please advise on the technical aspect of the question. Maple easily says that the answer is (1/2,1/2) but it is not clear why there is such an answer.

@sand15 

Thanks a million! I try to follow the steps in your file and understand how it is calculated. Thanks again.

@mmcdara 

Thank you so much! This was very helpful. But I could not apply it to my system of equations because it is not polynomial. I put my model. You see that the answer is obtained as (1/2, 1/2). But we don't know how Maple obtained it or whether there is an alternate solution. 

I uploaded my Maple code.

restart;

Z := (beta__c*y__c^2 + beta__m*y__c*y__n + beta__n*y__n^2)/(2*h*(-r__c*y__c - r__n*y__n + h)) + (beta__c*(1 - y__c)^2 + beta__m*(1 - y__c)*(1 - y__n) + beta__n*(1 - y__n)^2)/(2*h*(h - (1 - y__c)*r__c - (1 - y__n)*r__n));
sys := {diff(Z, y__c) = 0, diff(Z, y__n) = 0};
solve(sys, [y__c, y__n]);

                     [[y__c =1/2-, y__n =1/2]]

Thank you so much. It was my mistake. I had missed the second term in the function. I corrected it. The second term is similar to the first one but all x changes to 1-x. Thank you!

@acer

Thank you so much. This is the output of lprint

w__1 = 3/2*beta*((1+beta)^2-1/alpha)/u(1-alpha)^2;
w__2 = (1+beta)^2/u;
w__3 = (psol-c)/(psol-s);
Delta = -4*w__1*w__3+w__2^2;
p = -1/2*(-2*alpha*u*v*w__1+ Delta^(1/2)*s-Delta^(1/2)*v-w__2*s+w__2*v)/w__1/alpha/u;

@Thomas Richard 

Thanks million. It is very useful indeed.

@acer 

Thank you! I tried to copy-paste the code, but its appearance was very unclear. I will simplify it and submit to main part:

w1 := (3*beta*((1 + beta)^2 - 1/alpha))/(2*u(1 - alpha)^2);
w2 := (1 + beta)^2/u + ((beta^2/(u*(1 - alpha))) . ((-delta*v + psol - sExp)/(psol - s)));
w3 := (psol - c + beta*(-delta*v + psol - sExp))/(psol - s);
Delta := -4*w1*w3 + w2^2;

p := max(solve((w2 - sqrt(Delta))/(2*w1) = alpha*u*(v - psol)/(v - s), psol));

@Thomas Richard 

Thank you so much. Would you please show me an example to see how to use 'fnormal' with 'solve'? Shall I call any package through 'with'? 

@Mac Dude 

Thank you so much. I solve a function as 'solve(f,x)=0'. There are multiple roots and I need to consider the largest root. The issue is that I need to work with a maximum of multiple roots. So, I use this answer 'max(solve(f,x)=0)'. Then, if the maximum of roots is complex and the software considers it blank, that's perfect for me. The issue is that when applying the max function, I receive this error: "Error, complex argument to max/min"