I am having an issue with the functions that I have defined using the solutions given by NDSolve. The easiest way to explain my issue is through the code itself, so here it is:
boundary1 = 0.3;
boundary2 = 1;
h[r_] = 10^(-2);
\[Mu]1 = 8*10^4;
\[Mu]2 = 8*10^4;
p = 1300;
nsolution1[(w_)?NumericQ] := 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, 1}][[1]];
radialstretch[(w_)?NumericQ, (r_)?NumericQ] := lr[r] /. nsolution1[w];
thetastretch[(w_)?NumericQ, (r_)?NumericQ] := lt[r] /. nsolution1[w]
angle[(w_)?NumericQ, (r_)?NumericQ] := \[Beta][r] /. nsolution1[w]
radialstretchderivative[(w_)?NumericQ, (r_)?NumericQ] := D[lr[u] /. nsolution1[w], u] /. u -> r
thetastretchderivative[(w_)?NumericQ, (r_)?NumericQ] := D[lt[u] /. nsolution1[w], u] /. u -> r
R[(w_)?NumericQ, (r_)?NumericQ] := r*thetastretch[w, r]
Rderivative[(w_)?NumericQ, (r_)?NumericQ] := thetastretch[w, r] + r*thetastretchderivative[w, r]
f1[(w_)?NumericQ, (r_)?NumericQ] := -(radialstretch[w, r]^2 - Rderivative[w, r]^2)^(1/2)
Z[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[Re[f1[w, t]], {t, 0.001, r}]
When I try to evaluate any of the functions I've defined at specific values of r, I do get an answer for them, however, when I try plotting the functions, that does not work. I understand it has to do with the way Mathematica is evaluating everything, however, I do not understand it well enough to know where exactly the problem is. Here is an example:
Plot[Evaluate[lr[r] /. nsolution1[1.2]], {r, 0.001, 1}]
This gives me a plot that is correct, but if I do the following:
Plot[radialstretch[1.2, r], {r, 0.001, 1}]
which should give me the same exact plot, instead I get the error: "NDSolve::dsvar: 0.0010204081428571428` cannot be used as a variable". together with a bunch of other ones.
Also, I thought either ParametricNDSolve or ParametricNDSolveValue might be better options in my case since I have w as a variable, however, that is somehow giving me a different answer for the functions when compared to the answers using NDSolve. I might have to make a separate question for that because I am extremely confused as to why I would get a different answer.
