# Why code is running and not giving any output for NSolve?

I am trying to solve the given equation in the code and want to plot 'N1' with 'P1'. But the code is not giving me any output and still showing running. What is the problem with it? If anyone can figure it out will be appreciated.

dela = -1.5;(*Subscript[\[CapitalDelta], L] ; Laser-Lower polariton \
detuning*)
g1 = 0.0003; (* Subscript[g, 0] ; vacuum optomechanical strength *)
del0 = -0.006; (* \[Delta] ;  exciton-photon detuning*)
ome = 3.014; (* \[CapitalOmega] ; Rabi splitting *)
omem = 0.125; (* Subscript[\[CapitalOmega], m] ; mechanical \
resonator's frequency *)
kexc = 0.002;  (* Subscript[\[Kappa], exc  ;  ]exciton decay rate*)
kk = .2;   (* \[Kappa]; cavity decay rate *)
gma = 0.00001;  (* \[CapitalGamma] ; phonon dcay rate *)
nth = 0;
kal = (kk + kexc)/2 - del0*(kk - kexc)/(2*ome);
g0 = g1*(1 - del0/ome)/2;
kb = 0.1;
eq1 = -I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*(Conjugate[a1] - a1) -
I*(P1*g1/(2*Sqrt[2]*ome))*(alab - Conjugate[alab] +
Conjugate[betab] - betab) + kal*N1 +
kb*(alab*a1 + Conjugate[alab]*Conjugate[a1] +
betab*Conjugate[a1] +
Conjugate[betab]*
a1 + (N1 - 2*(Abs[a1])^2)*(b1 + Conjugate[b1])) == 0;
eq2 = I*dela*alab + I*omem*alab -
I*g0*((2*N1 - Abs[a1]^2 - Np + 2*2*(Abs[b1])^2 + b1^2)*
Conjugate[a1] - (Conjugate[betab] + 2*alab)*b1 -
alab*Conjugate[b1]) + I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 -
I*(P1*g1/(2*Sqrt[2]*ome))*(N1 - Np + (Conjugate[a1])^2 -
b1^2) + (kal + gma)*alab/2 +
kb/2*(Conjugate[betab]*b1 +
alab*Conjugate[b1] + (Np - 2*(Abs[b1])^2)*Conjugate[a1] +
2*alab*b1 - b1^2*Conjugate[a1]) == 0;
eq3 = -I*g0*((N1 - 2*(Abs[a1])^2)*(Conjugate[b1] - b1) +
Conjugate[alab]*Conjugate[a1] - alab*a1 + Conjugate[betab]*a1 -
betab*Conjugate[a1]) -
I*(P1*g1/(2*Sqrt[2]*ome))*(Conjugate[betab] - alab +
Conjugate[alab] - betab) + gma*Np - gma*nth == 0;
eq4 = -I*dela*betab + I*omem*betab -
I*g0*((2*N1 - (Abs[a1])^2 + Np - 2*(Abs[b1])^2)*
a1 - (2*betab - a1*b1 + Conjugate[alab])*b1 +
betab*Conjugate[b1]) - I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 -
I*(P1*g1/(2*Sqrt[2]*ome))*(N1 + Np + a1^2 + b1^2 + 1) + (kal +
gma)*betab/2 +
kb/2*((2*betab - a1*b1 + Conjugate[alab])*
b1 + (Np - 2*(Abs[b1])^2)*a1 + betab*Conjugate[b1]) == 0;
eq5 = -I*dela*a1 - I*g0*(Conjugate[alab] + betab) -
I*(P1/Sqrt[2]*(1 - del0/(2*ome))) -
I*(P1*g1/(2*Sqrt[2]*ome))*(b1 + Conjugate[b1]) + kal*a1/2 +
kb*(Conjugate[alab] + betab)/2 == 0;
eq6 = I*omem*b1 - I*g0*N1 -
I*(P1*g1/(2*Sqrt[2]*ome))*(a1 + Conjugate[a1]) + gma*b1/2 == 0;
sol[P1_?NumericQ] := FindRoot[{eq1, eq2, eq3, eq4, eq5, eq6}, {N1, alab, Np, betab, a1,
b1}];
Plot[{Evaluate[N1] /. sol[P1]}, {P1, 0, 20}]

• You have not really read and taken into account any comments made to you in this question why-this-code-ndsolve-is-not-working And it will really help if you write the ODE's each in separate equation, then use these later inside NDSolve instead of writing everything all in one command. No need to do that. Makes it hard for you and others to see things. The same here. Better write each equation on its own (say eq1, eq2 etc..).then use these inside the Solve commands to make it more clear. May 19, 2020 at 4:04
• @Nasser thanks for making corrections.
– vini
May 19, 2020 at 4:06
• At a minimum you need to assign a numeric value to P1. Are there any known constraints on {N1, alab, Np, betab, a1, b1}? May 19, 2020 at 5:19
• Edit your question to show what you are doing. We shouldn't have to guess. It appears that you used sol without including its argument. N1 /. sol[P1] May 19, 2020 at 5:49
• ...initialization...;tbl=Table[{P1,N1}/.NMinimize[Abs[..]+Abs[..]+Abs[..]+Abs[..]+Abs[..]+Abs[..],{N1,alab,Np,betab,a1,b1}][[2]],{P1,0,20,1}];ListPlot[tbl,Joined->True] pretty quickly shows a plot of P1 versus N1 for the range of P1 with no warnings or errors when I substitute in the left hand side of your six equations inside each of those Abs[..] and make no other changes. Adjust the step size of P1 as desired. Looking at the full result from the NMinimize can show you how well it approximated zero.
– Bill
May 19, 2020 at 7:18

I had thought that the skeleton of the code I showed above in my comment would be enough for you to be able to fill in the details. Apparently that wasn't true. So here is all the code.

dela=-1.5;g1=0.0003;del0=-0.006;ome=3.014;omem=0.125;kexc=0.002;kk=.2;gma=0.00001;nth=0;
kal=(kk+kexc)/2-del0*(kk-kexc)/(2*ome);g0=g1*(1-del0/ome)/2;kb=0.1;

expr1=-I*(P1/Sqrt[2]*(1-del0/(2*ome)))*(Conjugate[a1]-a1)-I*(P1*g1/(2*Sqrt[2]*ome))*(alab-
Conjugate[alab]+Conjugate[betab]-betab)+kal*N1+kb*(alab*a1+Conjugate[alab]*Conjugate[a1]+
betab*Conjugate[a1]+Conjugate[betab]*a1+(N1-2*(Abs[a1])^2)*(b1+Conjugate[b1]));

expr2=I*dela*alab+I*omem*alab-I*g0*((2*N1-Abs[a1]^2-Np+2*2*(Abs[b1])^2+b1^2)*Conjugate[a1]-
(Conjugate[betab]+2*alab)*b1-alab*Conjugate[b1])+I*(P1/Sqrt[2]*(1-del0/(2*ome)))*b1-I*(P1*
g1/(2*Sqrt[2]*ome))*(N1-Np+(Conjugate[a1])^2-b1^2)+(kal+gma)*alab/2+kb/2*(Conjugate[betab]*
b1+alab*Conjugate[b1]+(Np-2*(Abs[b1])^2)*Conjugate[a1]+2*alab*b1-b1^2*Conjugate[a1]);

expr3=-I*g0*((N1-2*(Abs[a1])^2)*(Conjugate[b1]-b1)+Conjugate[alab]*Conjugate[a1]-alab*a1+
Conjugate[betab]*a1-betab*Conjugate[a1])-I*(P1*g1/(2*Sqrt[2]*ome))*(Conjugate[betab]-
alab+Conjugate[alab]-betab)+gma*Np-gma*nth;

expr4=-I*dela*betab+I*omem*betab-I*g0*((2*N1-(Abs[a1])^2+Np-2*(Abs[b1])^2)*a1-(2*betab-a1*
b1+Conjugate[alab])*b1+betab*Conjugate[b1])-I*(P1/Sqrt[2]*(1-del0/(2*ome)))*b1-I*(P1*g1/(2*
Sqrt[2]*ome))*(N1+Np+a1^2+b1^2+1)+(kal+gma)*betab/2+kb/2*((2*betab-a1*b1+Conjugate[alab])*
b1+(Np-2*(Abs[b1])^2)*a1+betab*Conjugate[b1]);

expr5=-I*dela*a1-I*g0*(Conjugate[alab]+betab)-I*(P1/Sqrt[2]*(1-del0/(2*ome)))-I*(P1*g1/(2*
Sqrt[2]*ome))*(b1+Conjugate[b1])+kal*a1/2+kb*(Conjugate[alab]+betab)/2;

expr6=I*omem*b1-I*g0*N1-I*(P1*g1/(2*Sqrt[2]*ome))*(a1+Conjugate[a1])+gma*b1/2;

tbl=Table[{P1,N1}/.NMinimize[Abs[expr1]+Abs[expr2]+Abs[expr3]+Abs[expr4]+Abs[expr5]+Abs[expr6],
{N1,alab,Np,betab,a1,b1}][[2]],{P1,0,20,1}];

ListPlot[tbl,Joined->True]


Now, unless I've made a mistake, I have made no change to any of your six equations other than removing the right hand side so that I can NMinimize the sum of the Abs of each of those.

To be very clear, trying to do NMinimize[Abs[expression]] is not the same as Solve[expression==0]. The first one of those tries to get close to expression==0.

So, as I wrote above, you should also look carefully at

Table[NMinimize[Abs[expr1]+Abs[expr2]+Abs[expr3]+Abs[expr4]+Abs[expr5]+Abs[expr6],
{N1,alab,Np,betab,a1,b1}],{P1,0,20,1}]


and see how well it did. The first item in each of those 21 results will be the minimum of the sum that it was able to find. You might like all those first items to be exactly zero, but they won't be. The question is: Are the sums good enough to give you an idea what your plot looks like.

There is another issue. Plot might do several hundred or even several thousand iterations, finding the value of the N1 at each of those points. This only does 21 points. So it can easily be that there is behavior between any pair of those points which is not displayed in the resulting ListPlot. Changing the step size for P1 might help you get a little better idea what the resulting graph will be like, but without taking hundreds or a thousand times longer to display a result.

There is another issue. You almost certainly intended that alab,Np,betab,a1,b1 could have complex values, not just real values. If you read the documentation for NMinimize and you click on the orange Details and Options you will find: By default, all variables are assumed to be real. Perhaps you can explore changing that default or looking at other functions that might let you find better minima if you search among complex values for alab,Np,betab,a1,b1 and still be able to rapidly get an estimate of what your graph looks like.

All this was ONLY intended to rapidly give you some idea what your graph might look like.

Hopefully this is enough to show you how to use this cheap reasonably fast trick that I intended.

If you can quickly determine whether this is good enough for your purposes then that was the only goal. If it isn't good enough then perhaps letting Solve run for a few hours or days or months might be enough to get you the original Plot that you desired.

Rationalize all of your constants so that you can specify WorkingPrecision later to use arbitrary-precision.

dela = -3/2;(*Subscript[Δ,L];Laser-Lower polariton detuning*)
g1 = 3*^-4;(*Subscript[g,0];vacuum optomechanical strength*)
del0 = -6*^-3;(*δ;exciton-photon detuning*)
ome = 3014*^-3;(*Ω;Rabi splitting*)
omem = 1/8;(*Subscript[Ω,m];mechanical resonator's frequency*)
kexc = 1/500;(*Subscript[κ,exc;]exciton decay rate*)
kk = 1/5;(*κ;cavity decay rate*)
gma = 10^-5;(*Γ;phonon dcay rate*)
nth = 0;
kal = (kk + kexc)/2 - del0*(kk - kexc)/(2*ome);
g0 = g1*(1 - del0/ome)/2;
kb = 1/10;
eq1 = -I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*
(Conjugate[a1] - a1) -
I*(P1*g1/(2*Sqrt[2]*ome))*
(alab - Conjugate[alab] + Conjugate[betab] - betab) +
kal*N1 + kb*(alab*a1 + Conjugate[alab]*Conjugate[a1] +
betab*Conjugate[a1] + Conjugate[betab]*a1 +
(N1 - 2*(Abs[a1])^2)*(b1 + Conjugate[b1])) == 0;
eq2 = I*dela*alab + I*omem*alab -
I*g0*((2*N1 - Abs[a1]^2 - Np + 2*2*(Abs[b1])^2 + b1^2)*
Conjugate[a1] - (Conjugate[betab] + 2*alab)*b1 - alab*Conjugate[b1]) +
I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 -
I*(P1*g1/(2*Sqrt[2]*ome))*(N1 - Np + (Conjugate[a1])^2 - b1^2) + (kal +
gma)*alab/2 +
kb/2*(Conjugate[betab]*b1 +
alab*Conjugate[b1] + (Np - 2*(Abs[b1])^2)*Conjugate[a1] + 2*alab*b1 -
b1^2*Conjugate[a1]) == 0;
eq3 = -I*g0*((N1 - 2*(Abs[a1])^2)*(Conjugate[b1] - b1) +
Conjugate[alab]*Conjugate[a1] - alab*a1 + Conjugate[betab]*a1 -
betab*Conjugate[a1]) -
I*(P1*g1/(2*Sqrt[2]*ome))*(Conjugate[betab] - alab + Conjugate[alab] -
betab) + gma*Np - gma*nth == 0;
eq4 = -I*dela*betab + I*omem*betab -
I*g0*((2*N1 - (Abs[a1])^2 + Np - 2*(Abs[b1])^2)*
a1 - (2*betab - a1*b1 + Conjugate[alab])*b1 + betab*Conjugate[b1]) -
I*(P1/Sqrt[2]*(1 - del0/(2*ome)))*b1 -
I*(P1*g1/(2*Sqrt[2]*ome))*(N1 + Np + a1^2 + b1^2 + 1) + (kal + gma)*
betab/2 +
kb/2*((2*betab - a1*b1 + Conjugate[alab])*b1 + (Np - 2*(Abs[b1])^2)*a1 +
betab*Conjugate[b1]) == 0;
eq5 = -I*dela*a1 - I*g0*(Conjugate[alab] + betab) -
I*(P1/Sqrt[2]*(1 - del0/(2*ome))) -
I*(P1*g1/(2*Sqrt[2]*ome))*(b1 + Conjugate[b1]) + kal*a1/2 +
kb*(Conjugate[alab] + betab)/2 == 0;
eq6 = I*omem*b1 - I*g0*N1 - I*(P1*g1/(2*Sqrt[2]*ome))*(a1 + Conjugate[a1]) +
gma*b1/2 == 0;

eqns = {eq1, eq2, eq3, eq4, eq5, eq6};


Using a modified version of what Bill suggested for NMinimize along with using arbitrary-precision (specified WorkingPrecision).

data = Table[{N1, P1} /.
NMinimize[Total[Abs[First[#]]^2 & /@ eqns], {N1, alab, Np, betab, a1, b1},
WorkingPrecision -> 15][[2]], {P1, 0, 20, 1/2}];

ListLinePlot[data]


Note that in your revised code you used the wrong syntax for FindRoot. See the documentation. You must supply initial estimates for each variable.