Solving set of nonlinear equations with multiple roots

I would like to solve the set of nonlinear equations for x & y, where the set of equations if defined in the function "Uniaxial", given below:

X=fsolve(@(X) Uniaxial(X,paramz),X0);
function F = Uniaxial(X,paramz)
    J0=paramz(1); V=paramz(2); P=paramz(3);
    x=X(1); y=X(2); x2=x^2; y2=y^2; x4=x2^2; y4=y2^2;       
    F(1) = (V^2)-(1/x2/y2)*(J0*(x4*y2-1)/(x2*y2*(J0-x2-y2+3)-1));
    F(2) = (V^2)-(1/x2/y2)*(J0*(x2*y4-1)/(x2*y2*(J0-x2-y2+3)-1)-y*P);  
end

V and P are two parameters that change in the ranges [0,1] and [0,2], respectively. J0 is a constant that is set equal to 100. I try to solve the problem by running the values of parameters V and P in a loop and finding the value of y, which is one of the unknowns of the set. If the solution were achived properly, I would be able to polot the value of V against y (the seconf unknown) as in the figure below (Fig.6 in Extreme Mechanics Letters 38 (2020) 100752):

The blue and black curves represetns solution for two different values of P.

Here, it can be seen that whereas for certain values of V there are three solution (black and blue points) for some other values of V slightly larger than 0.5, no solution exists, which becomes evident when I run my program: exactly at the first maximum of the function V(y), the solutions do not match this curve.

Natrually, one way to go around this problem is to change the initial values x0 after this point, but the main problem here is that I don't know of any way to know this point a priori.

Is there a way to solve this problem?