# Is there any other way to find roots of a non-monotonic function?

I have the relation of $k$ and $w$, and now I need to plot $k$ vs $w$($k\sim 10^7$, $w\sim 10^{13}$). The relation is given as follows:

parameters:

\[Epsilon][2] = \[Epsilon][1] = 9.5; \[Epsilon][3] = 1;\[Mu] = 1.257*10^-6;c = 3.*10^8;\[Epsilon]0 = 8.85*10^-12;d = 20.*10^-9;
\[Sigma] = Im[1/(-I* w)*Abs[(7.5*10^14*(1.6*10^-19)^2)/(0.2*0.91*10^-30)]];


relation:

1 - ((\[Epsilon][1] Sqrt[-\[Epsilon][2] w^2/c^2 + k^2] - \[Epsilon][
2] Sqrt[-\[Epsilon][1] w^2/c^2 +
k^2] - \[Sigma] Sqrt[-\[Epsilon][1] w^2/c^2 + k^2]
Sqrt[-\[Epsilon][2] w^2/c^2 +
k^2]/(\[Epsilon]0  w ))/(\[Epsilon][
1] Sqrt[-\[Epsilon][2] w^2/c^2 + k^2] + \[Epsilon][
2] Sqrt[-\[Epsilon][1] w^2/c^2 +
k^2] - \[Sigma] Sqrt[-\[Epsilon][1] w^2/c^2 + k^2]
Sqrt[-\[Epsilon][2] w^2/c^2 +
k^2]/(\[Epsilon]0  w ))) (\[Epsilon][
3] Sqrt[-\[Epsilon][2] w^2/c^2 + k^2] - \[Epsilon][
2] Sqrt[-\[Epsilon][3] w^2/c^2 + k^2])/(\[Epsilon][
3] Sqrt[-\[Epsilon][2] w^2/c^2 + k^2] + \[Epsilon][
2] Sqrt[-\[Epsilon][3] w^2/c^2 + k^2])Exp[-2 d Sqrt[-\[Epsilon][2] w^2/c^2 + k^2]]=0


The method ContourPlot doesn't work well, so I turn to use FindRoot and NSolve to find the list of solutions $(k,w)$ so I can plot it manually. However, even FindRoot doesn't work well. For example, when I set w=3 Pi * 10 ^12, the solution of k should be $k=1.41179*10^8$. However, I cannot use FindRoot to find this solution since the LHS of the relation is a non-monotonic function: Now I have to choose the value of k0 very carefully for FindRoot to start with, or Mathematica cannot give me the correct answer. However, what I need to do is to plot the relation between $k$ and $w$, so I really don't want to find k0 by looking into the plot with different $w$ one by one, I hope to find a way to make Mathematica smartly choose k0 and find the solution.

I am grateful for any help, thank you!

The standard way is to enlist NDSolve to solve your equation as a DAE:

ksol = First@
NDSolveValue[{
obj == 0 /. k -> k[w],   (* obj = OP's "relation" w/o "= 0" *)
k[3 Pi*10^12] == (k /. FindRoot[obj /. w -> 3 Pi*10^12, {k, 1.4*^8}]),
y'[w] == 0, y[3 Pi*10^12] == 0},      (* dummy ODE *)
{k}, {w, 1*^12, 3 Pi*10^12}, StartingStepSize -> 1];
(* N.B. step size >> 3 Pi*10^12 * \$MachineEpsilon *)

ListLinePlot[ksol, PlotRange -> All]


If greater precision is desired, you can use ksol to estimate initial values for FindRoot.

• excellent, thank you very much! – Jieyu You Nov 22 '17 at 4:01
• Is there any essential difference between NDSolve and NDSolveValue? – Jieyu You Nov 22 '17 at 4:09
• @JieyuYou NDSolveValue[ODE, y, {x, a, b}, options] and y /. NDSolve[ODE, y, {x, a, b}, options] are completely equivalent. The difference is the form of the return value and how you use it. I tend to prefer NDSolve, esp. for systems with multiple dependent variable or multiple solutions. More and more I see experienced users saying they prefer NDSolveValue, which is convenient for DEs of the form y'[x] == f[x, y[x]]. – Michael E2 Nov 22 '17 at 12:54