Neural Network for a highly nonlinear relation (not a function)

Hello and greetings everybody.

I intend to model a neural network (MLP preferably), for an input-output relation as the following figure which is not a function and is highly nonlinear

I tried a lot of strategies (Neuron size increase/decrease, layer change, different methods like RBF and...) all of which lead to absolute failure. This figure in physics is called "A bistable relation" . I know that by some rotation, I can invert this figure to functional relation and then it would be a piece of cake to model yet I am strongly prohibited to do this. How can I model a neural network with an accpetable MSE for such relation. I already attached the file for data creation and the default number of samples is k = 10000. It would be kind of you to help me at this. Thanks in advance for your time devoted to this equation.

%%%Part I ==> Constants Input
format long 

g2 = 3.5155;                  
g3 = g2;                     
g1 = 0;                                     
CP = 200;                      
Oc1 =  0;                   
k = 10000;                        %Iteration time


%%Part II => Bistability Relation, Im and Re part of Rho21 based on
% different values of Oc (D21=D32=0)

             
%%Preallocating Matrices

G21 = zeros(1,k);

G2 = zeros(1,k);

G3 = zeros(1,k);

G1 = zeros(1,k);

G31 = zeros(1,k);

G32 = zeros(1,k);

OC1 = zeros(1,k);

OC2 = zeros(1,k);

OC3 = zeros(1,k);

rho1 = zeros(1,k);

rho2 = zeros(1,k);

rho3 = zeros(1,k);

x = zeros(1,k);

C = zeros(1,k);


wp = zeros(1,k);
  %Detuning range (typically for Rho21)
    
    Y = linspace(0,50,k);       %Outputs
  


for j = 1 :k
    
    %Substituting varibales into preallocating matrices
    
    G2(1,j) = g2;
    
    G3(1,j) = g3;
    
    G1(1,j) = g1;
    
    G21(1,j) = (G1(1,j) + G2(1,j)) ./ 2;
    
    G32(1,j) = (G3(1,j) + G2(1,j)) ./ 2;
    
    G31(1,j) = (G3(1,j) + G1(1,j)) ./ 2;
    
    OC1(1,j) = Oc1;
    
    
    C(1,j) = CP;
    
    %Rho21 Derivation for various Oc's
    
    Ro1 = [-G2(1,j) -(1i*Y(1,j)) 0 (1i*Y(1,j)) 0 0 0 -G2(1,j);
        -(2i*Y(1,j)) (-1i*wp(1,j)-G21(1,j)) -(1i*OC1(1,j)) 0 0 0 0 -(1i*Y(1,j));
        0 -(1i*OC1(1,j)) (-1i*wp(1,j)-G31(1,j)) 0 (1i*Y(1,j)) 0 0 0;
        (2i*Y(1,j)) 0 0 (1i*wp(1,j)-G21(1,j)) 0 (1i*OC1(1,j)) 0 (1i*Y(1,j));
        (1i*OC1(1,j)) 0 (1i*Y(1,j)) 0 -G32(1,j) 0 0 (2i*OC1(1,j));
        0 0 0 (1i*OC1(1,j)) 0 (1i*wp(1,j)-G31(1,j)) -(1i*Y(1,j)) 0;
        -(1i*OC1(1,j)) 0 0 0 0 -(1i*Y(1,j)) -G32(1,j) -(2i*OC1(1,j));
        0 0 0 0 (1i*OC1(1,j)) 0 -(1i*OC1(1,j)) -G3(1,j)];
    
    B1 = [-G2(1,j);-(1i*Y(1,j));0;(1i*Y(1,j));(1i*OC1(1,j));0;-(1i*OC1(1,j));0];
   
    
    R1 = Ro1 \ B1;
    
    rho1(1,j) = R1(4);
    
    %input-output relation : |x| in terms of |y|
 
   x(1,j) = (2 .* Y(1,j)) - (1i .* C(1,j) .* rho1(1,j));
    

end

X = abs(x);         %Inputs


%%%Part III ==> Plotting

figure
plot(X,Y)
xlabel('input |y|')
ylabel('output |x|') 
hold on
grid on