# Obtaining different answers when using NDSolve vs ParametricNDSolveValue

I am obtaining two different answers for a curve when solving the same exact system of differential equations when using NDSolve vs when using ParametricNDSolveValue. Here it is when using NDSolve:

h[r_] = 10^(-2);
boundary1 = 0.5;
boundary2 = 1;
\[Mu]1 = 8*10^4;
\[Mu]2 = 8*10^4;
p = 1300;
w = 1.48323;
nsolution1 = NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h[r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) -
\[Mu]1*h[r]*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h[r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) +
\[Mu]1*h[r]*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w,
lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r, 0.001, boundary1}][[1]];
thetastretch[(r_)?NumericQ] := lt[r] /. nsolution1
angle[(r_)?NumericQ] := \[Beta][r] /. nsolution1
radialstretchderivative[(r_)?NumericQ] := D[lr[u] /. nsolution1, u] /. u -> r
thetastretchderivative[(r_)?NumericQ] := D[lt[u] /. nsolution1, u] /. u -> r
R[(r_)?NumericQ] := r*thetastretch[r]
Rderivative[(r_)?NumericQ] := thetastretch[r] + r*thetastretchderivative[r]
Z[(r_)?NumericQ] := NIntegrate[Re[f1[t]], {t, 0.001, r}]
p1 := ParametricPlot[{R[r], Z[r] + 1}, {r, 0.001, boundary1}, PlotStyle -> Green]
nsolution2 = Quiet[NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]2*h[r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) -
\[Mu]2*h[r]*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]2*h[r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) +
\[Mu]2*h[r]*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]2*h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0,
lt[boundary1] == thetastretch[boundary1], lr[boundary1] == g[boundary1], \[Beta][boundary1] == angle[boundary1]}, {lr, lt, \[Beta]}, {r, boundary1, boundary2}][[1]]];
thetastretch2[(r_)?NumericQ] := Quiet[lt[r] /. nsolution2]
angle2[(r_)?NumericQ] := Quiet[\[Beta][r] /. nsolution2]
radialstretchderivative2[(r_)?NumericQ] := Quiet[D[lr[u] /. nsolution2, u] /. u -> r]
thetastretchderivative2[(r_)?NumericQ] := Quiet[D[lt[u] /. nsolution2, u] /. u -> r]
R2[(r_)?NumericQ] := r*thetastretch2[r]
Rderivative2[(r_)?NumericQ] := thetastretch2[r] + r*thetastretchderivative2[r]
Z2[(r_)?NumericQ] := Quiet[NIntegrate[Re[f2[t]], {t, boundary1, r}]]
sol2 = Solve[Z2[boundary1] + x2 == Z[boundary1], Plus[x2]];
x2 = x2 /. sol2[[1]];
y2 := 1 + x2;
p2 := ParametricPlot[{R2[r], Z2[r] + y2}, {r, boundary1, boundary2}, PlotStyle -> Red]
Show[p1, p2, PlotRange -> All]


Now, here it is solving it using ParametricNDSolveValue:

h = 1/10^2;
\[Mu]1 = 8*10^4;
\[Mu]2 = 8*10^4;
p = 1300;
boundary1 = 0.5;
boundary2 = 1;
nsolution1 = ParametricNDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) -
\[Mu]1*h*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0,
(-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]1*h*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) +
(2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*h*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w, lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r, 0.001, boundary1}, {w}, AccuracyGoal -> 10, PrecisionGoal -> 10];
radialstretch1[(w_)?NumericQ, (r_)?NumericQ] := lr[w][r] /. nsolution1
thetastretch1[(w_)?NumericQ, (r_)?NumericQ] := lt[w][r] /. nsolution1
angle1[(w_)?NumericQ, (r_)?NumericQ] := \[Beta][w][r] /. nsolution1
radialstretchderivative1[(w_)?NumericQ, (r_)?NumericQ] := D[lr[w][t] /. nsolution1, t] /. t -> r
thetastretchderivative1[(w_)?NumericQ, (r_)?NumericQ] := D[lt[w][t] /. nsolution1, t] /. t -> r
R1[(w_)?NumericQ, (r_)?NumericQ] := r*thetastretch1[w, r]
Rderivative1[(w_)?NumericQ, (r_)?NumericQ] := thetastretch1[w, r] + r*thetastretchderivative1[w, r]
f1[(w_)?NumericQ, (r_)?NumericQ] := -Sqrt[radialstretch1[w, r]^2 - Rderivative1[w, r]^2]
Z1[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[f1[w, r], {t, 0.001, r}]
nsolution2 = ParametricNDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]2*h*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) -
\[Mu]2*h*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0,
(-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]2*h*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]2*h*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) +
(2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]2*h*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[boundary1] == thetastretch1[w, boundary1], lr[boundary1] == g1[w, boundary1], \[Beta][boundary1] == angle1[w, boundary1]}, {lr, lt, \[Beta]}, {r, boundary1, boundary2}, {w}, AccuracyGoal -> 10,
PrecisionGoal -> 10];
radialstretch2[(w_)?NumericQ, (r_)?NumericQ] := lr[w][r] /. nsolution2
thetastretch2[(w_)?NumericQ, (r_)?NumericQ] := lt[w][r] /. nsolution2
angle2[(w_)?NumericQ, (r_)?NumericQ] := \[Beta][w][r] /. nsolution2
radialstretchderivative2[(w_)?NumericQ, (r_)?NumericQ] := D[lr[w][t] /. nsolution2, t] /. t -> r
thetastretchderivative2[(w_)?NumericQ, (r_)?NumericQ] := D[lt[w][t] /. nsolution2, t] /. t -> r
R2[(w_)?NumericQ, (r_)?NumericQ] := r*thetastretch2[w, r]
Rderivative2[(w_)?NumericQ, (r_)?NumericQ] := thetastretch2[w, r] + r*thetastretchderivative2[w, r]
f2[(w_)?NumericQ, (r_)?NumericQ] := -Sqrt[radialstretch2[w, r]^2 - Rderivative2[w, r]^2]
Z2[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[f2[w, r], {t, boundary1, r}]
p1 := ParametricPlot[{R1[1.48323, r], Z1[1.48323, r] + 1}, {r, 0.001, boundary1}, PlotStyle -> Green]
Clear[x2];
Clear[y2];
sol2 = Solve[Z2[1.48323, boundary1] + x2 == Z1[1.48323, boundary1], Plus[x2]];
x2 = x2 /. sol2[[1]];
y2 = 1 + x2;
p2 := ParametricPlot[{R2[1.48323, r], Z2[1.48323, r] + y2}, {r, boundary1, boundary2}, PlotStyle -> Red]
Show[p1, p2, PlotRange -> All]


As you can see, the equations are exactly the same, the parameters used are exactly the same, however, the results are different. I would love to use ParametricNDSolveValue because it is much more convenient for my code, which is much longer than this, and it works much faster, but I need to first know why it is giving me a different answer. If anyone has any insight, please let me know.

• There are many typos in the code, for example, h[r] is not defined for nsolution1. Commented Jan 29, 2023 at 5:59
• It was in my notebook, I somehow missed copying it here. I edited it now. Commented Jan 29, 2023 at 7:03

Solutions computed with NDSolve and ParametricNDSolve are the same as it follows from the test

Clear["*"]

boundary1 = 0.5;
boundary2 = 1;
\[Mu]1 = 8*10^4;
\[Mu]2 = 8*10^4;
p = 1300;
h = 1/10^2;

nsolutionp1 =
ParametricNDSolve[{lt[r] + r*Derivative[1][lt][r] ==
lr[r]*Cos[\[Beta][r]],
p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h*
Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*
Derivative[1][\[Beta]][r]*
Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*
Sin[\[Beta][
r]]*(Derivative[1][lr][
r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) ==
0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h*
Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]1*h*r*
Derivative[1][\[Beta]][r]*
Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*
Cos[\[Beta][
r]]*(Derivative[1][lr][
r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*
h*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w,
lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r,
0.001, boundary1}, {w}];
lr[w][r] /. nsolutionp1;
thetastretch1[(w_)?NumericQ, (r_)?NumericQ] := lt[w][r] /. nsolutionp1;
angle1[(w_)?NumericQ, (r_)?NumericQ] := \[Beta][w][r] /. nsolutionp1;
D[lr[w][t], t] /. t -> r /. nsolutionp1;
thetastretchderivative1[(w_)?NumericQ, (r_)?NumericQ] :=
D[lt[w][t], t] /. t -> r /. nsolutionp1;
R1[(w_)?NumericQ, (r_)?NumericQ] := r*thetastretch1[w, r];
Rderivative1[(w_)?NumericQ, (r_)?NumericQ] :=
thetastretch1[w, r] + r*thetastretchderivative1[w, r];
f1[(w_)?NumericQ, (r_)?NumericQ] := -Sqrt[
Z1[(w_)?NumericQ, (r_)?NumericQ] :=
NIntegrate[f1[w, r], {t, 0.001, r}];
g1[(w_)?NumericQ, (r_)?
NumericQ] := (1/
thetastretch1[w, r]^2)) + ((\[Mu]1 -
2); nsolutionp2 =
ParametricNDSolve[{lt[r] + r*Derivative[1][lt][r] ==
lr[r]*Cos[\[Beta][r]],
p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]2*h*
Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*
Derivative[1][\[Beta]][r]*
Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*
Sin[\[Beta][
r]]*(Derivative[1][lr][
r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) ==
0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]2*h*
Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]2*h*r*
Derivative[1][\[Beta]][r]*
Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*
Cos[\[Beta][
r]]*(Derivative[1][lr][
r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]2*
h*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0,
lt[boundary1] == thetastretch1[w, boundary1],
lr[boundary1] == g1[w, boundary1], \[Beta][boundary1] ==
angle1[w, boundary1]}, {lr, lt, \[Beta]}, {r, boundary1,
boundary2}, {w}];


To test solutions we plot

{test1 =
Plot[Evaluate[{lr[w][r], lt[w][r], \[Beta][w][r]} /.
w -> 1.48323 /. nsolutionp1], {r, 0.001, boundary1},
PlotLegends -> {lr, lt, \[Beta]}],
test2 = Plot[
Evaluate[{lr[w][r], lt[w][r], \[Beta][w][r]} /. w -> 1.48323 /.
nsolutionp2], {r, boundary1, boundary2},
PlotLegends -> {lr, lt, \[Beta]}]}


NDSolve code

Clear[{h, w, f1}]

boundary1 = 0.5;
boundary2 = 1;
\[Mu]1 = 8*10^4;
\[Mu]2 = 8*10^4;
p = 1300;
w = 1.48323; h[r_] := 1/10^2;

nsolution1 =
NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]],
p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h[r]*
Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*
Derivative[1][\[Beta]][r]*
Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*
Sin[\[Beta][
r]]*(Derivative[1][lr][
r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) ==
0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h[r]*
Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]1*h[r]*r*
Derivative[1][\[Beta]][r]*
Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*
Cos[\[Beta][
r]]*(Derivative[1][lr][
r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*
h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w,
lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r,
0.001, boundary1}][[1]];
thetastretch[(r_)?NumericQ] := lt[r] /. nsolution1
angle[(r_)?NumericQ] := \[Beta][r] /. nsolution1
D[lr[u] /. nsolution1, u] /. u -> r
thetastretchderivative[(r_)?NumericQ] :=
D[lt[u] /. nsolution1, u] /. u -> r
R[(r_)?NumericQ] := r*thetastretch[r]
Rderivative[(r_)?NumericQ] :=
thetastretch[r] + r*thetastretchderivative[r]
Z[(r_)?NumericQ] := NIntegrate[Re[f1[t]], {t, 0.001, r}]
p1 := ParametricPlot[{R[r], Z[r] + 1}, {r, 0.001, boundary1},
PlotStyle -> Green]
g[r_] := (1/
thetastretch[r]^2)) + ((\[Mu]1 -
nsolution2 =
NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]],
p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]2*h[r]*
Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*
Derivative[1][\[Beta]][r]*
Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*
Sin[\[Beta][
r]]*(Derivative[1][lr][
r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) ==
0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]2*h[r]*
Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]2*h[r]*r*
Derivative[1][\[Beta]][r]*
Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*
Cos[\[Beta][
r]]*(Derivative[1][lr][
r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]2*
h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0,
lt[boundary1] == thetastretch[boundary1],
lr[boundary1] == g[boundary1], \[Beta][boundary1] ==
angle[boundary1]}, {lr, lt, \[Beta]}, {r, boundary1,
boundary2}][[1]];
thetastretch2[(r_)?NumericQ] := Quiet[lt[r] /. nsolution2]
angle2[(r_)?NumericQ] := Quiet[\[Beta][r] /. nsolution2]
Quiet[D[lr[u] /. nsolution2, u] /. u -> r]
thetastretchderivative2[(r_)?NumericQ] :=
Quiet[D[lt[u] /. nsolution2, u] /. u -> r]
R2[(r_)?NumericQ] := r*thetastretch2[r]
Rderivative2[(r_)?NumericQ] :=
thetastretch2[r] + r*thetastretchderivative2[r]
f2[(r_)?NumericQ] :=
Z2[(r_)?NumericQ] := Quiet[NIntegrate[Re[f2[t]], {t, boundary1, r}]]
sol2 = Solve[Z2[boundary1] + x2 == Z[boundary1], x2];
x2 = x2 /. sol2[[1]];
y2 := 1 + x2;
p2 := ParametricPlot[{R2[r], Z2[r] + y2}, {r, boundary1, boundary2},
PlotStyle -> Red]


Now we can compare two solutions as follows

{Show[test1,
Plot[Evaluate[{lr[r], lt[r], \[Beta][r]} /. nsolution1], {r, 0.001,
boundary1}, PlotLegends -> {lr, lt, \[Beta]},
PlotStyle -> {{Black, Dashed}, {Black, Dashed}, {Black,
Dashed}}]],
Show[test2,
Plot[Evaluate[{lr[r], lt[r], \[Beta][r]} /. nsolution2], {r,
boundary1, boundary2}, PlotLegends -> {lr, lt, \[Beta]},
PlotStyle -> {{Black, Dashed}, {Black, Dashed}, {Black, Dashed}}]]}


Therefore, all discrepancies coming from definitions p1,p2 and other functions. For example, Z and Z1 are differ, while definitions f1 are not differ

We should look what NIntegrate doing in two cases. Finally we found out typos in definitions Z1[w, r], Z2[w,r]

Z1[(w_)?NumericQ, (r_)?NumericQ] :=
NIntegrate[f1[w, r], {t, 0.001, r}];
Z2[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[f2[w, r], {t, boundary1, r}]


Obviously it should be

Z1[(w_)?NumericQ, (r_)?NumericQ] :=
NIntegrate[f1[w, t], {t, 0.001, r}];
Z2[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[f2[w, t], {t, boundary1, r}]


After correction we have same picture for NDSolve and ParametricNDSolve (dashed lines)

• Thank you for this. Can you see some obvious mistake that I have made when defining the functions using the solutions from ParametricNDSolve? I don't know another way of defining them for them to work as when I do so from NDSolve. Commented Jan 29, 2023 at 19:41
• More precisely, x2 in two definitions is not the same since Z1[1.48323, boundary1]=-0.431163, while Z[boundary1]=-0.230469. Commented Jan 30, 2023 at 2:06
• Thanks for the additional help. Since f1 are not different, all the issues have to be due to the integrating of it. Somehow in the second case(the ParametricNDSolve one), using NIntegrate on f1 is doing something different. I will try some stuff and post an update soon. Commented Jan 30, 2023 at 5:34
• @juv95 There is a typo in definition Z1[w, r] since you used under NIntegrate function f1[w,r] instead of f1[w,t]`. Commented Jan 30, 2023 at 6:19
• Ahh, that was so careless of me... thank you for pointing it out. Everything is working out well now, except I had to make a change on the definition of Z1. Initially, I had it defined as NIntegrate[f1[w,t],{t,0.001,r}], and that was not working. When I switched it to NIntegrate[f1[w,r],{r,0.001,t}] instead, it works great. I am not sure why that is the case, but I'll mention it here in case it helps anyone in the future. Commented Jan 30, 2023 at 7:34