# Evaluation functions with lower and upper limit ODE?

How do I find the numerical value of y at 0 and also at 1. Here is my code. I appreciate any help.

    ϵ = 10^-4;
R = Rationalize[1.5472, 0];
sol = ParametricNDSolveValue[{y''[r] + 3 y'[r]/r == Sinh[y[r]],
y[ϵ] == y0, y'[ϵ] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + 32001/40]}, {y, y'},
{r, ϵ, R}, {y0}, Method -> "StiffnessSwitching", WorkingPrecision -> 30];
ff = FindRoot[Last[sol[y0]][R], {y0, -12.25, -11.}, Evaluated -> False][[1, 2]]


## 1 Answer

You can avoid the singularity problems near r==0! Change the boundarie conditions to y[R]==0,y'[R]==ys

R = Rationalize[1.5472, 0];
Clear[ys];
Y = ParametricNDSolveValue[{r y''[r] + 3 y'[r] == r Sinh[y[r]],
y[R] == 0, y'[R] == ys,WhenEvent[r == 1, y'[r] -> y'[r] + 32001/40]}, y, {r, 0, R},{ys} ,Method -> "StiffnessSwitching"]


and force the solution to fullfill Y'[ys]==0

ys=ys/. NMinimize[ Y'[ys]^2, ys][]
(* 0.674335*)


The functionvalues evaluate to

{Y[ys] ,Y[ys]}
(*{-2008.16, -0.764241}*)


at r==0 and r==1

Hopefully that's what you're looking for!