How to plot multicolor phase daigram with three parameters?


enter image description here

I am trying to plot like the diagram given above. I have two nonlinear equations that I have separated in imaginary and complex forms (u1,v1, and x1,y1). As shown in the figure in the x-axis I want ‘p0’ which varies from (0,2.5*10^3) and in the y-axis, I want ‘del0’ which varies from (0,3) and inside the plot where rectangular color region box it is ‘NL=u1^2+v1^2’which varies from (0,6*10^7). I am trying this but unable to plot this type of diagram and gives me an error: “solve was unable to solve the system with inexact coefficients. The \ the answer was obtained by solving a corresponding exact system and \ numericizing the result. >>” Should I give another command to get this type of plot or anything else? If anyone can short it out will be appreciated.

A1 = 1; B1 = 0; delta = -1.5;(*Subscript[Δ, L] ; Laser-Lower polariton \ detuning*) g0 = .315; (* Subscript[g, 0] ; vacuum optomechanical strength *) ome = 3.014; (* Ω ; Rabi splitting *) ome1 = .125; (* Subscript[Ω, m] ; mechanical \ resonator's frequency *) kexc = 0.002;  (* Subscript[κ, exc  ;  ]exciton decay rate*) kk = .2;   (* κ; cavity decay rate *) gma = 0.00001;  (* Γ ; phonon dcay rate *) ka = (kk + kexc)/2 - del0*(kk - kexc)/(2*ome); Sol = NSolve[{-delta*v1 - ka*u1/2 - (kk - kexc)*B1*g0*x1*u1/ome == 0,      delta*u1 + g0*(1 - del0/ome)*A1*x1*u1 +        p0*(1 - del0/(2*ome))/Sqrt[2] + p0*g0*B1*x1/(Sqrt[2]*ome) -        ka*v1/2 - (kk - kexc)*g0*x1*B1*v1/ome == 0,      ome1*y1 - gma*x1/2 ==       0, -ome1*x1 + g0/2*(1 - del0/ome)*A1*(u1^2 + v1^2) +        p0*g0*u1*B1/(Sqrt[2]*ome) - gma*y1/2 == 0}, {u1, v1, x1, y1},     del0]; ContourPlot[{Evaluate[(u1^2 + v1^2)] /. Sol}, {p0, 0, 5}, {del0, 0,    3}, PlotLegends -> BarLegend[{"LakeColors", {0, 6}}]]  After running this code I am getting the plot. Which is not the same as above. I don't know what the problem with it?