Is there any way to join two ODE solutions so that I can use the final product as a function like an usual ODE solution?
My problem is that I have a function α(ϕ) defined in two regimes, so I have to solve two ODEs. The first starts with physical initial conditions at ϕ=0, the second starts when the first stops and continues until ϕ=2. The first step is then to join the two ODE solutions from ϕ=[0,2], the second to symmetrise this function so that, in the end, the solution is valid for ϕ=[-2,2].
Below is a MWE in which I solve the ODEs, then extract the solution and time in arrays and use append! to glue the solutions from the two regimes, and finally create a DiffEqArray to relate the complete solution with the values of ϕ. Plotting this array show the correct solution, but I cannot use it as a function, which is necessary for me to solve another ODE. Is there an alternative to what I have tried to do? Maybe reconstructing an ODESolution type from the two ODE solutions?
using DifferentialEquations, Plots, RecursiveArrayTools
const global lp=1.
s2 = 81.0
d = -0.0009
B = 1.0
β = 1.0
k = 1.0
α₀ = 0.083163
ϕ₀ = 2.0
function quantum!(dv,v,p,ϕ)
s2,d,B,β,k=p
α = v[1]
#@show (ϕ*s2*sin(2*d*v[1])+2*d*sinh(s2*v[1]*ϕ))/(-v[1]*s2*sin(2*d*v[1])+2*d*cos(2*d*v[1])+2*d*cosh(s2*v[1]*ϕ))
dv[1] = (ϕ*s2*sin(2*d*α)+2*d*sinh(s2*α*ϕ))/(-α*s2*sin(2*d*α)+2*d*cos(2*d*α)+2*d*cosh(s2*α*ϕ))
end
# When dα/dϕ = 1, we reach the classical regime, and stop the integration
condition(v,ϕ,integrator) = (ϕ*s2*sin(2*d*v[1])+2*d*sinh(s2*v[1]*ϕ))/(-v[1]*s2*sin(2*d*v[1])+2*d*cos(2*d*v[1])+2*d*cosh(s2*v[1]*ϕ))==1.0
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition,affect!)
# Initial Condition at the bounce
a₀ = [α₀]
classspan = (0,ϕ₀)
probquant = ODEProblem(quantum!,a₀,classspan,(s2,d,B,β,k))
solquant = solve(probquant,Tsit5(),callback=cb)
function classic!(du,u,p,ϕ)
αc = u[1]
dαc = u[2]
du[1] = dαc
du[2] = 3*(-dαc^3/sqrt(2)+dαc^2+dαc/sqrt(2)-1)
end
# IC retrieved from end of quantum phase
init = [last(solquant);1.0]
classspan = (last(solquant.t),ϕ₀)
probclass = ODEProblem(classic!,init,classspan)
solclass = solve(probclass,Tsit5())
# I tried to join the two regime like this:
# α(ϕ) for ϕ>0
solu = append!(solquant[1,:],solclass[1,:]) # α
solt = append!(solquant.t,solclass.t) # ϕ
# α(ϕ) for ϕ<0
soloppu = reverse(solu)
soloppt = -reverse(solt)
pop!(soloppu) # Remove -0.0 (redundant)
pop!(soloppt)
# Join the two solutions
soltotu = append!(soloppu,solu)
soltott = append!(soloppt,solt)
soltot = DiffEqArray(soltotu,soltott)
plot(soltot) # Plot is fine, but how to use soltot as a function?