# Joining two ODE solutions
**URL:** https://discourse.julialang.org/t/joining-two-ode-solutions/63710
**Category:** General Usage
**Tags:** diffeq, ode
**Created:** [June 28, 2021, 6:06pm UTC](https://discourse.julialang.org/t/joining-two-ode-solutions/63710 "2021-06-28T18:06:15Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@Free\_ion](https://discourse.julialang.org/u/Free_ion)
#### Post date: [June 28, 2021, 6:06pm UTC](https://discourse.julialang.org/t/joining-two-ode-solutions/63710/1 "2021-06-28T18:06:15Z")
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?
```julia
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?
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 2, 2021, 12:13pm UTC](https://discourse.julialang.org/t/joining-two-ode-solutions/63710/2 "2021-07-02T12:13:04Z")
@Vaibhavdixit02 didn’t you have a solution for this?
---
### Author:  [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)
#### Post date: [July 2, 2021, 2:02pm UTC](https://discourse.julialang.org/t/joining-two-ode-solutions/63710/3 "2021-07-02T14:02:27Z")
Yeah it’s this code in DiffEqParameEstim.jl [DiffEqParamEstim.jl/multiple\_shooting\_objective.jl at 227ececc2af8b24746e395c88973ad9f7f87e0af · SciML/DiffEqParamEstim.jl · GitHub](https://github.com/SciML/DiffEqParamEstim.jl/blob/227ececc2af8b24746e395c88973ad9f7f87e0af/src/multiple_shooting_objective.jl#L58-L64)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 2, 2021, 3:08pm UTC](https://discourse.julialang.org/t/joining-two-ode-solutions/63710/4 "2021-07-02T15:08:40Z")
We should get something on this added to SciMLBase.
---
### Author:  [@Free\_ion](https://discourse.julialang.org/u/Free_ion)
#### Post date: [August 10, 2021, 1:30pm UTC](https://discourse.julialang.org/t/joining-two-ode-solutions/63710/5 "2021-08-10T13:30:59Z")
Sorry for the late reply!
So it seems that I should add a line similar to
```julia
sol_new = DiffEqBase.build_solution(prob, alg, sol_loss.t, sol_loss.u, retcode =:Success)
```
from Vaibhavdixit02’s answer. I’m a bit confused on the `prob` option though. Applying this to my case, I initially thought I could have used something like
```julia
soltot = DiffEqBase.build_solution(prob, Tsit5(), soltott, soltotu,retcode = :Success)
```
instead of the `DiffEqArray` at the end of my code. However, I don’t know what I should use instead of prob since I have both `probquant` and `probclass` as problems. This issue does not seem to appear in Vaibhavdixit02’s answer since `prob` is passed as an argument in `function multiple_shooting_objective`. Any idea?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 26, 2021, 4:31am UTC](https://discourse.julialang.org/t/joining-two-ode-solutions/63710/6 "2021-08-26T04:31:07Z")
The prob is used for the `f` and `p` in the interpolation of the lazy algorithms. For the most part, you should be okay merging like this.