Possible bug: Wrong differentiation of ParametricFunction inside Parallelize?

Question

In the note RootSearch problem when called from within ParallelDo I described the problem with Ted Ersek's RootSearch package and came to the conclusion that it is caused by the fact that the derivative with respect to the parameter from the ParametricFunction, which is obtained as a solution to the differential equation using the ParametricNDSolve function, is incorrectly calculated inside the ParallelDo loop.

The following example illustrates the problem. First, we solve an ordinary differential equation and formulate the result of the solution in the form of a function ϕ1p, which depends on the z coordinate and the parameters β and Λ. We also define the ϕ1p1b function as the derivative of the ϕ1p function with respect to the β parameter. This function and its numerical value at β=0.4 will be compared with the ϕ1p derivative with respect to the β parameter, which we will calculate inside the Do and ParallelDo loops. This derivative within the cycles is denoted as LHS1b.

Inside the Do loop, the ϕ1p1b and LHS1b functions and their numerical values are the same. However, inside the ParallelDo loop, the resulting expression for the LHS1b function is gigantic, so the Short form is used to print it. Mathematica cannot calculate the numerical value of LHS1b at all.

Λs = {2, 5,
  10}; (* list of values for ParallelDo cycle *)
Km = 8; k = \[Infinity]; \[Mu] = 1; \[Nu] = 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 = { (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, 1000}}];
ϕ1p[z_, β_, Λ_] = ϕ[β, \
Λ][z] /. sol1p;
Print["ϕ1p[z,β,Λ]=", ϕ1p[
   z, β, Λ]];
ϕ1p1z[z_, β_, Λ_] =
 D[ϕ1p[z, β, Λ], z]
ϕ1p1b[z_, β_, Λ_] =
 D[ϕ1p[z, β, Λ], β]
(*=====================================================================*)


DistributeDefinitions[ϕ1p, ϕ1p1z, ϕ1p1b, BC1];
βm = 0.4;
"Do cycle ===================== "
Do[
  LHS[β_] = ϕ1p[1, β, Λm];
  Print["LHS[β]=", LHS[β],
   "\nϕ1p[1,β,Λm]=", ϕ1p[
    1, β, Λm]];

  LHS1b[β_] =
   D[ϕ1p[1, β, Λm], β];
  Print["LHS1b[β]=", Short[LHS1b[β], 3],
   "\nϕ1p1b[1,β,Λm]=", ϕ1p1b[
    1, β, Λm]];
  Print["LHS1b[βm]=", Short[LHS1b[βm], 3],
   "\nϕ1p1b[1,βm,Λm]=", ϕ1p1b[
    1, βm, Λm]];

  RHS[β_] = BC1[β, Km, Λm];
  Print["D[RHS[β],β]=", D[RHS[β], β]];
  , {Λm, Λs}
  ];

"ParallelDo cycle ============= "
ParallelDo[
  LHS[β_] = ϕ1p[1, β, Λm];
  Print["LHS[β]=", LHS[β],
   "\nϕ1p[1,β,Λm]=", ϕ1p[
    1, β, Λm]];

  LHS1b[β_] =
   D[ϕ1p[1, β, Λm], β];
  Print["LHS1b[β]=", Short[LHS1b[β], 3],
   "\nϕ1p1b[1,β,Λm]=", ϕ1p1b[
    1, β, Λm]];
  Print["LHS1b[βm]=", Short[LHS1b[βm], 3],
   "\nϕ1p1b[1,βm,Λm]=", ϕ1p1b[
    1, βm, Λm]];

  RHS[β_] = BC1[β, Km, Λm];
  Print["D[RHS[β],β]=", D[RHS[β], β]];
  , {Λm, Λs}
  ];

It seems to me that the described problem is a bug in the Mathematica system, although I admit that I am doing something wrong. If so, where am I wrong?

EDIT #1: Same error occurs when computing derivative over coordinate z.