RootSearch problem when called from within ParallelDo

Question

I have been using Ted Ersek's RootSearch package for Wolfram Mathematica for many years. The package worked great, including inside the Do loop. However, inside the ParallelDo loop, it issues a diagnostic message, for example:

Thread :: tdlen: Objects of unequal length in {0,0} + {0,0} + <<8>> + <<9>> cannot be combined.

At the same time, the calculation result seems to remain correct.

At the moment I am not ready to post minimal working example here as my code is quite huge (Currently its size is 1k). But I noted that Ted Ersek is active at this forum so I hope that he let me know how I could send him my code.

Here is an MWE derived from my code. I am guessing that the problem has something to do with the fact that the equation for RootSearch uses the ParametricNDSolve solution.

Needs["Ersek`RootSearch`"];

debug = 0; (* controls how extensive is debug printout *)

Λs = {1.1, 2, 5, 10, 50, 100, 500}; (* list of values for ParallelDo cycle *)

results = {{"k", "μ", "ν", "K", "Λ", "βcrit1", "βcrit2", "$KernelID"}};
SetSharedVariable[results];
DistributeDefinitions[RootSearch];

Km = 8; k = ∞; μ = 1; ν = 2; (* particular parameters from outer Do cycle not shown here *)

"Boundary condition at z=1: BC1[β_,K_,Λ_]" (* used in RootSearch *)
BC1[β_, K_, Λ_] = -(((-1 + K) (3 K^2 - 
     2 β) β ϕ[1])/(
  2 (K^2 - β) (K^2 - β + K^2 Λ)))

"Equation to be solved: eqn="
eqn = {1/256 (49 β (1 - (
         64 (8 - 7 z)^2 β (128 - (8 - 7 z)^2 β))/(512 - 
           8 (8 - 7 z)^2 β)^2) ϕ[z] + 
      1/(512 - 
         8 (8 - 7 z)^2 β)^2 6272 β (32 (-1 + 
            z) z^3 β^2 - z^4 β^2 + 
         64 z^2 β (1 - 6 (-1 + z)^2 β) + 
         1024 (-1 + z) z β (-1 + 2 (-1 + z)^2 β) + 
         4096 (-2 + (-1 + z)^2 β - (-1 + 
               z)^4 β^2)) ϕ[z] + 
      56 (8 - 7 z) β Derivative[1][ϕ][z] + 
      4 (64 - (8 - 7 z)^2 β + 
         64 Λ) (ϕ′′[z]) == 
   0, ϕ[0] == 1, 
  Derivative[1][ϕ][0] == (7 β (-3 + 2 β) ϕ[0])/(
   16 (-1 + β) (-1 + β - Λ))}

"Solution of ParametricNDSolve:"
sol1p = ParametricNDSolve[
  eqn, ϕ, {z, 0, 1}, {{β, 0, 1}, {Λ, 1, 10000}}]
ϕ1p[z_, β_, Λ_] = ϕ[β, Λ][z] /. sol1p
ϕ1p1[z_, β_, Λ_] = D[ϕ1p[z, β, Λ], z]

(*===================================================================*)

DistributeDefinitions[ϕ1p, ϕ1p1, BC1];
ParallelDo[
  ClearAll[βcrit1, βcrit2, Err1, Err2];
  If[debug > 2
   , Print@
     Plot[{ϕ1p1[1, β, Λm], 
       BC1[β, 
         Km, Λm] /. {ϕ[1] -> ϕ1p[
           1, β, Λm]}}, {β, 0, 0.9999}
      , AxesLabel -> {"β", {ϕ', BC1}}
      , PlotLegends -> {ϕ', BC1}
      , PlotLabel -> 
       StringForm[
        "k=``, μ=``, ν=``, K=``, Λ=``", 
        k, μ, ν, Km, Λm]];
   ];
  
  βcrit1[$KernelID] = 
   Quiet@FindRoot[ϕ1p1[1, β, Λm] == 
       BC1[β, 
        Km, Λm] /. {ϕ[1] -> ϕ1p[
         1, β, Λm]}, {β, 0.999, 0.1, 0.9999}
     ];
  
  Err1[$KernelID] = (ϕ1p1[1, β, Λm] - 
       BC1[β, Km, Λm]) /. {ϕ[1] -> ϕ1p[
        1, β, Λm]} /. βcrit1[$KernelID];
  If[debug > 0, 
   Print@ToString@
     StringForm[
      "FindRoot1:: k=``, μ=``, ν=``, K=``, \
Λ=``, β1=``, Err1=``", k, μ, ν, 
      Km, Λm, Quiet@βcrit1[$KernelID][[1, 2]], 
      Err1[$KernelID]]];
  
  Print@StringForm[
    "6, kernel ``/``, Err1=``", $KernelID, $KernelCount, 
    Err1[$KernelID]];
  
  βcrit2[$KernelID] = 
   RootSearch[ϕ1p1[1, β, Λm] == 
      BC1[β, 
       Km, Λm] /. {ϕ[1] -> ϕ1p[
        1, β, Λm]}, {β, 0.1, 0.9999}
    , InitialSamples -> 10
    , RootTest -> (Abs[#2] < 10^-4 &)
    ];
  Print@StringForm[
    "7, kernel ``/``, βcrit2=``", $KernelID, $KernelCount, \
βcrit2[$KernelID]];
  If[Length[βcrit2[$KernelID]] > 0
   , Err2[$KernelID] = (ϕ1p1[1, β, Λm] - 
        BC1[β, Km, Λm]) /. {ϕ[
         1] -> ϕ1p[1, β, Λm]} /. 
     First[βcrit2[$KernelID]]
   , Err2[$KernelID] =.;
   ];
  
  Print@StringForm[
    "8, kernel ``/``, Err2=``", $KernelID, $KernelCount, 
    Err2[$KernelID]];
  If[debug > 0, 
   Print@ToString@
     StringForm[
      "RootSearch2:: k=``, μ=``, ν=``, β2=``, Err2=``", 
      k, μ, ν, Quiet@βcrit2[$KernelID][[1, 1, 2]], 
      Err2[$KernelID]
      ]
   ];
  
  If[Not@NumericQ[Err1[$KernelID]] || 
    Abs[Err1[$KernelID]] > 
     10^-4, βcrit1[$KernelID] = {β -> β1}];
  If[Not@NumericQ[Err2[$KernelID]] || 
    Abs[Err2[$KernelID]] > 
     10^-4, βcrit2[$KernelID] = {{β -> β2}}];
  AppendTo[results, 
   elem = {k, μ, ν, 
     Km, Λm, βcrit1[$KernelID][[1, 
       2]], βcrit2[$KernelID][[1, 1, 2]], $KernelID}];
  Print["elem=", elem];
  
  , {Λm, Λs}
  ];
(*===================================================================*)

results // TableForm
results

EDIT #2: I have located a chunck of code in ResearchRoot package that yields the warning message. It is located inside a defintion of RootSearch function. In the code below the source of warning is located between two Print commands added by me:

Print["myRootSearch 5: x=", x];
If[invalidQ,$Failed,
  (*else*)f=Evaluate[{SamplePoints[#],lhs-rhs/.x->#}]&;
  f=ReplacePart[f,Set,{1,0}];
  Print["myRootSearch 5.1: f=",f];
  u=Function[Evaluate[
        withTemp[{setTemp[y,(lhs-rhs)/.x->#],
            setTemp[yyp,((lhs-rhs)/D[lhs-rhs,x])/.x->#]},
          Hold[
            SamplePoints[#]=y,
            Which[TrueQ[y==0],0,NumericQ[yyp],yyp
              ,True,Indeterminate]
            ]]]];
  Print["myRootSearch, 5.2"];
  u=ReplacePart[
      u/.{withTemp->With,setTemp->Set},{CompoundExpression},{{1,
          2,0}},{{1}}];

These Prints produce the following output (from Kernel 2):

myRootSearch 5: x=[Beta]

myRootSearch 5.1: f=(KiamyRootSearchPrivate`SamplePoints[#1]=(ParametricFunction[Expression: [Phi] Parameters: {#1,500}

][#1,500]^[Prime])[1]+(7 (192-2 #1) #1 ParametricFunction[Expression: [Phi] Parameters: {#1,500}

][#1,500][1])/(2 (64-#1) (32064-#1))

Thread::tdlen : Objects of unequal length in {0,0}+{0,0}+<<8>>+<<9>> cannot be combined.

myRootSearch, 5.2

As you see, the warning goes between print points 5.1 and 5.2.

EDIT #3:

I have found that error is caused by D[lhs-rhs,x].