I have a function that reads: f[P_?NumericQ, Ob_?NumericQ] := -Sqrt[P/Ob^2 + Sqrt[P^2 + Ob^2 w^2]/Ob^2] Tan[ Sqrt[-(P/Ob^2) + Sqrt[P^2 + Ob^2 w^2]/Ob^2]] + Sqrt[-(P/Ob^2) + Sqrt[P^2 + Ob^2 w^2]/Ob^2] Tanh[Sqrt[P/Ob^2 + Sqrt[P^2 + Ob^2 w^2]/Ob^2]] For each pair of numerical values of P and Ob, there are infinitely many solutions for w (denoted as w0, w1, w2, w3, ...). I want to write a concise code that does the following: 1. Loops over values of Ob from 0 to 1. 2. For each value of Ob, finds a critical value of P (Pcr) that makes the ratio of w2/w1 approximately equal to 3 (with reasonable tolerance). 3. Stores or outputs these critical values of P for each Ob. I tried but couldn't come up with a successful one. Thanks for any help in advance. Below is one of my trials: Remove["Global`*"]; tolerance = 0.05; results = {}; For[ObValue = 0.01, ObValue <= 1, ObValue += 0.01, Pcr = {}; roots = NSolve[f[P, ObValue] == 0 && 0 < w < 20, w, Reals]; If[Length[roots] > 1, w1 = w /. roots[[1]]; w2 = w /. roots[[2]]; Pcr = FindRoot[Abs[w2/w1 - 3] <= tolerance, {P, .1}]; If[NumberQ[P /. Pcr], AppendTo[results, {ObValue, P /. Pcr}];];];] results