Hi again, some of my relevant code and images, to illustrate my question above:

```
BesselJInitAmount = 1.75
BesselYInitAmount = 0.37
FnOrder = [1.27]
tspanForward = (5.0, 20.0)
tspanBackward = (5.0, 0.01)
function OrdBesselDiffEqToSolve!(ddu,du,u,p,t)
OrderParam = p[1]
ddu[1] = -(du[1]/t) - (1 - ((OrderParam/t)^2))*u[1]
end
function OrdBesselInitCondCalc(BesJamt, BesYamt, Ord, tVal)
StartVal = (BesJamt*besselj(Ord,tVal)) + (BesYamt*bessely(Ord,tVal))
StartDeriv = ((BesJamt/2)*(besselj((Ord-1),tVal) - besselj((Ord+1),tVal))) +
((BesYamt/2)*(bessely((Ord-1),tVal) - bessely((Ord+1),tVal)))
return [StartVal, StartDeriv]
end
function OrdBesselTheorySoln(t, BesJamt, BesYamt, Ord)::BigFloat
return ((BesJamt*besselj(Ord,t)) + (BesYamt*bessely(Ord,t)))
end
InitialConditions = OrdBesselInitCondCalc(BesselJInitAmount, BesselYInitAmount,
FnOrder[1], tspanForward[1])
xInit = [InitialConditions[1]]
vInit = [InitialConditions[2]]
probFor = SecondOrderODEProblem(OrdBesselDiffEqToSolve!,vInit,xInit,tspanForward,FnOrder)
probBack = SecondOrderODEProblem(OrdBesselDiffEqToSolve!,vInit,xInit,tspanBackward,FnOrder)
#
#
#@time solFor = solve(probFor, dense=true, reltol=1e-12, abstol=1e-12)
#@time solBack = solve(probBack, dense=true, reltol=1e-12, abstol=1e-12)
#
#@time solFor = solve(probFor, Tsit5(), dense=true, reltol=1e-12, abstol=1e-12)
#@time solBack = solve(probBack, Tsit5(), dense=true, reltol=1e-12, abstol=1e-12)
#
@time solFor = solve(probFor, RKO65(), dense=true, reltol=1e-12, abstol=1e-12, dt=0.001)
@time solBack = solve(probBack, RKO65(), dense=true, reltol=1e-12, abstol=1e-12, dt=-0.001)
TvalsFor = getindex.(solFor.t, 1)
FprimeOfTFor = getindex.(solFor.u, 1)
FofTFor = getindex.(solFor.u, 2)
#
TvalsBack = getindex.(solBack.t, 1)
FprimeOfTBack = getindex.(solBack.u, 1)
FofTBack = getindex.(solBack.u, 2)
#
u_analyticalFor = @. OrdBesselTheorySoln(TvalsFor,
BesselJInitAmount, BesselYInitAmount, FnOrder[1])
errFor = abs.(u_analyticalFor - FofTFor)
@. errFor[errFor == 0 ] = NaN # set zeros to NaN to avoid problems caused by log(0) = -Inf
#
#
u_analyticalBack = @. OrdBesselTheorySoln(TvalsBack,
BesselJInitAmount, BesselYInitAmount, FnOrder[1])
errBack = abs.(u_analyticalBack - FofTBack)
@. errBack[errBack == 0 ] = NaN # set zeros to NaN to avoid problems caused by log(0) = -Inf
#
#
#
fineStepSize = 0.00299
fineForStepsWOends = BigInt(((tspanForward[2]-tspanForward[1]) ÷ fineStepSize)-1)
fineBackStepsWOends = BigInt(((tspanBackward[1]-tspanBackward[2]) ÷ fineStepSize)-1)
#
t_fineFor = append!(append!([tspanForward[1]],
[(tspanForward[1] + (stepN*fineStepSize)) for stepN = 1:fineForStepsWOends]),
[tspanForward[2]])
#
t_fineBack = append!(append!([tspanBackward[1]],
[(tspanForward[1] - (stepN*fineStepSize)) for stepN = 1:fineBackStepsWOends]),
[tspanBackward[2]])
#
fineForResids = [abs((solFor(time)[2] -
OrdBesselTheorySoln(time, BesselJInitAmount, BesselYInitAmount,
FnOrder[1]))) for time in t_fineFor]; @. fineForResids[fineForResids == 0] = NaN
#
fineBackResids = [abs((solBack(time)[2] -
OrdBesselTheorySoln(time, BesselJInitAmount, BesselYInitAmount,
FnOrder[1]))) for time in t_fineBack]; @. fineBackResids[fineBackResids == 0] = NaN
#
#
#
InterpDataPtResidPlots =
scatter(TvalsFor, errFor, <plot options>);
plot!(t_fineFor, fineForResids, <plot options>);
scatter!(TvalsBack, errBack, <plot options>);
plot!(t_fineBack, fineBackResids, <plot options>)
```

Hopefully this code snippet (and the attached image) clarifies my question…