How to get accurate solutions from FindRoot?

eq = 1 + θ/k^2 + (R (5 R + 4 R θ + R^3 θ - Sqrt[2] Sqrt[1/θ] (8 + 
4 (2 + R^2) θ + (2 + R^2)^2 θ^2) DawsonF[(R Sqrt[θ])/Sqrt[2]]))/(4 k^2);

n = 5;(*Increase the value to find more solution*)


FindInstance[eq == 0 && 1 < k < 3 && 0 < R < 1 && 1 < θ < 1000,
{R, k, θ}, n, RandomSeeding -> Automatic] // N // MatrixForm

enter image description here

For: n=15

n = 15;
FindInstance[
eq == 0 && 1 < k < 3 && 0 < R < 1 && 1 < θ < 1000, {R, 
k, θ}, n, RandomSeeding -> Automatic, 
WorkingPrecision -> 20] // MatrixForm

enter image description here

ContourPlot[Evaluate@Table[(eq /. θ -> j) == 0, {j, 1, 1000, 50}], {k, 1, 
3}, {R, 0, 1}, FrameLabel -> Automatic]

enter image description here

ContourPlot[Evaluate@Table[(eq /. θ -> j) == 0, {j, 1, 100}], {k, 1, 
3}, {R, 0, 1}, FrameLabel -> Automatic]

enter image description here


OP request: Solution by FindRoot (borrowed code form user: Henrik):

  Block[{R, k, \[Theta]}, 
  eq = {R, k, \[Theta]} \[Function] 
  Evaluate@
  Simplify[
  1 + \[Theta]/k^2 + 
   R^2/k^2 (1 + \[Theta] (1 + 1/(4 \[Theta]) + R^2/4) + (Sqrt[
         2/\[Theta]] 1/
         R (1 + \[Theta] (1 + (R^2)/2) + \[Theta]^2/
            2 (1 + (R^2)/2)^2) (-2 DawsonF[
           Sqrt[\[Theta]/2] R]))), \[Theta] > 0];]
     solr = {k, \[Theta]} \[Function] 
 Block[{R}, R /. FindRoot[eq[R, k, \[Theta]], {R, 1/10, 3/4}, Method -> "Secant"]] // Quiet;
 data = Block[{\[Theta] = 10}, Table[{k, solr[k, \[Theta]]}, {k, 1, 3, 1/10}]] // MatrixForm

enter image description here

As you can see values for R is not in domain: 0<R<1.


The division by R in the function eq is superfluous and makes problems in the root finding. Fortunately, some simplification can get you rid of that:

Block[{R, k, θ},
 eq = {R, k, θ} \[Function] Evaluate@Simplify[ 
      1 + θ/k^2 + R^2/k^2 (1 + θ (1 + 1/(4 θ) + R^2/4) + (Sqrt[ 2/θ] 1/ R (1 + θ (1 + (R^2)/2) + θ^2/ 2 (1 + (R^2)/2)^2) (-2 DawsonF[ Sqrt[θ/2] R]))),
      θ > 0
      ];
 ]

Second issue: For some parameter values of k, there may be more than one solution or no solutions at all. NSolve can find (hopefully) all solution:

solr = {k, θ} \[Function] Block[{R}, R /. NSolve[eq[R, k, θ] == 0 && 0 <= R <= 1, R, Reals, WorkingPrecision -> 25]]; 
data = Block[{θ = 100}, Table[{k, solr[k, θ]}, {k, 1, 3, 1/10}]]
ListPlot[Join @@ Thread /@ data, AxesLabel -> {"k", "R"}]

{{1, {0.1294234646803674672924767}}, {11/10, {0.1296843073043101973816057}}, {6/5, {0.1299709530763769620463354}}, {13/10, {0.1302836799112623098323127}}, {7/5, {0.1306227971663856528622202}}, {3/2, {0.1309886477076770723646643}}, {8/5, {0.1313816102517342127410522}}, {17/10, {0.1318021020175089797099938}}, {9/5, {0.1322505817265513300747067}}, {19/10, {0.1327275529978164798164777}}, {2,{0.1332335681913876276756329}}, {21/ 10, {0.1337692327655012292943036, 0.8928149869232539863495161}}, {11/5, {0.1343352102233942222258231, 0.7840644848720745156904206}}, {23/10, {0.1349322277412426656153432, 0.7061209842390900229330888}}, {12/5, {0.1355610825864964635669894, 0.6468555823332288287855214}}, {5/2, {0.1362226494580973891652682, 0.5999539682600664300172137}}, {13/5, {0.1369178889075194101436787, 0.5617652337478442737838912}}, {27/10, {0.1376478570337623102004883, 0.5300143374213614722559567}}, {14/5, {0.1384137166883060426333669, 0.5031957080178530201387233}}, {29/10, {0.1392167504801858849697920, 0.4802563992569981158981134}}, {3, {0.1400583759402635089617742, 0.4604211571053870525255899}}}

enter image description here