# FindRoot gives an imaginary value and it is outside Interpolated Function values

My code is

(*constants*)
c = 1;
G = 1;
R0 = 1.473;
k = 4.012/10^4;
γ = 2;
α = 0.2;
ϵ0 = (((1/k)*(R0/α)^γ)^(1 - γ))^(-((1 - γ)^2)^(-1));
ϵ0s = 1.603*10^38;
ScaleMeVfmtoMs = 8.92/10^7;
ScaleErgcmtoMs = 5.56/10^40;
β = ((4*Pi*ϵ0)/(k*ϵ0^(γ - 1))^(1/γ))*ScaleMeVfmtoMs;
(*Solving*)
TOV = Derivative[p][r] + (((α*M[r])/r^2)*p[r]^(1/γ)*(1 + p[r]/(p[r]/(k*ϵ0^(γ - 1)))^(1/γ))*(1 + (4*Pi*r^3*p[r])/(M[r]*c^2)))/
(1 - (2*G*M[r])/(c^2*r));
rTOV[p0_] := NDSolve[{TOV == 0, Derivative[M][r] == β*r^2*p[r]^(1/γ), p[0.001] == p0,
M[0.001] == (4/3)*Pi*0.001^3*((p[0.001]/(k*ϵ0^(γ - 1)))^(1/γ)/c^2)}, {p, M}, {r, 0.001, 30},
Method -> {"FixedStep", Method -> "ExplicitRungeKutta"}];
pTOV[r_, p0_] := p[r] /. rTOV[p0]
MTOV[r_, p0_] := M[r] /. rTOV[p0]
NR[p0_] := FindRoot[pTOV[r, p0] == 0, {r, 10}][[1,2]]
NM[p0_] := MTOV[NR[p0], p0][]

{NR[10^(-4)], NM[10^(-4)]}


.

The errors message are

FindRoot::nlnum: The function value {InterpolatingFunction[{{0.001,30.}},{5,15,1,{2013207},{4},0,0,0,0,Automatic,{},{},False},{{<<1>>}},{DeveloperPackedArrayForm,{<<1>>},{0.0001 +0. I,-1.4911*10^-8+0. I,0.0001 +0. I,-1.46601*10^-8+0. I,0.0001 +0. I,-1.44227*10^-8+0. I,0.0001 +0. I,<<37>>,0.0001 +0. I,-1.163*10^-8+0. I,0.0001 +0. I,-1.15577*10^-8+0. I,0.0001 +0. I,-1.14898*10^-8+0. I,<<4026364>>}},{Automatic}][<<1>>]} is not a list of numbers with dimensions {1} at {r} = {18.569 -6.99405*10^-7 I}.
FindRoot::nlnum: The function value {InterpolatingFunction[{{0.001,30.}},{5,15,1,{2013207},{4},0,0,0,0,Automatic,{},{},False},{{<<1>>}},{DeveloperPackedArrayForm,{<<1>>},{0.0001 +0. I,-1.4911*10^-8+0. I,0.0001 +0. I,-1.46601*10^-8+0. I,0.0001 +0. I,-1.44227*10^-8+0. I,0.0001 +0. I,<<37>>,0.0001 +0. I,-1.163*10^-8+0. I,0.0001 +0. I,-1.15577*10^-8+0. I,0.0001 +0. I,-1.14898*10^-8+0. I,<<4026364>>}},{Automatic}][<<1>>]} is not a list of numbers with dimensions {1} at {r} = {18.569 -6.99405*10^-7 I}.


This seems like some steps in FindRoot return a small imaginary value for next step and my function contains only real numbers. Is there any way to prevent FindRoot going to give this small number?

Also, the results are

{18.569 + 0. I, 0.903448 + 1.86409*10^-14 I}

• Your ODE has p[r]^(1/γ) where γ=2, so it's no surprise that complex numbers arise when you are trying to find when p[r] is zero. Probably using something like Abs[p[r]]^(1/γ) might help. – Carl Woll Jun 12 '19 at 16:43
• Hi Carl. I've edited as your suggestion. The code runs fine without any errors. But I'm not sure if doing this will affect the result? – Panithi Nakkhruea Jun 12 '19 at 16:56
• @PanithiNakkhruea You can use NR[p0_] := First[NSolve[pTOV[r, p0] == 0, r]][[1, 2]] // Quiet – Alex Trounev Jun 13 '19 at 4:59