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 #@show (ϕ*s2*sin(2*d*v)+2*d*sinh(s2*v*ϕ))/(-v*s2*sin(2*d*v)+2*d*cos(2*d*v)+2*d*cosh(s2*v*ϕ)) dv = (ϕ*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)+2*d*sinh(s2*v*ϕ))/(-v*s2*sin(2*d*v)+2*d*cos(2*d*v)+2*d*cosh(s2*v*ϕ))==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 dαc = u du = dαc du = 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?