# Solving chain of ODE in Julia

I am solving two different ODE whose solutions need to be matched. I am currently doing this by hand, which works great, but I would like to automatise this process. The second ODE takes one of its initial conditions and its initial time from the last values of the first solution. I put a MWE below. Now, I have naively tried to put together the first and second ODEs, as I show in the second code below (giving a StackOverflowError). It seems this is not the way to go, though. Does someone know how to do that?

EDIT: the second code is now the solution

using DifferentialEquations, Plots

# First ODE, which depends on some parameters

#ϕ₀ = abs(find_time!(s2,d))
#ϕ₀=2.0
s2,d=p
α = 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
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!)
α₀ = [a0]
tspan = (0,ϕ₀)
end

# Second ODE, which needs initial α and ϕ from the previous solution
function test!(dv,v,p,ϕ)
α = v[1]
dα = v[2]
dv[1] = dα
dv[2] = 3*(-dα^9/sqrt(2)+dα^8+dα^7/sqrt(2)-dα^6)
end
init = [1.06,1.0] # Put IC by hand
testspan = (0.75,3.0) # Put initial ϕ by hand
prob = ODEProblem(test!,init,testspan)
sol = solve(prob,Tsit5())
plot!(sol,vars=(0,1),xlims=(0,2))

using DifferentialEquations, Plots

s2,d=p
α = 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
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!)
α₀ = [a0]
tspan = (0,ϕ₀)

function classic!(du,u,p,ϕ)
αc = u[1]
dαc = u[2]
du[1] = dαc
du[2] = 3*(-dαc^9/sqrt(2)+dαc^8+dαc^7/sqrt(2)-dαc^6)
end
probclass = ODEProblem(classic!,init,classspan)
solclass = solve(probclass,Tsit5())

sol = DiffEqArray(solu,solt)
end
$$$$


last(soldadp) is an array, so [last(soldadp),1.0] is a Vector{Union{Vector{Float64},Float64}}, i.e. it's like [[2.0],1.0]. That's not what you wanted: you wanted [2.0,1.0] which means you wanted to use the concatenation operator ;. init = [last(soldadp);1.0] fixes your error.

Now, I don't know what you mean by sol = soldadp+solclass: those two time series are not aligned in time, and they hold different numbers of vectors. Did you want to append the plot of the second after the other? Then you can use plot! to mutate the plot. Or you can build an array that concatenates the two in the way you're interpreting, but the + operator doesn't seem to make much sense here. Either way, it's a 1 x N object and a 2 x M object, so it's going to need some user input on what you actually mean by "put these together".

• Thanks for pointing out the error! Indeed, sol as such does not make sense. What I really want as output is to extend the function α(ϕ) computed on [0,last(soldadp)] with the α(ϕ) computed on [last(soldadp), ϕ₀] (I will need this function later). In the end, I'm interested in u[1] only in the second ODE, so I thought I could just add the first component of the two solutions using sol = append!(soldadp[1,:],solclass[1,:]), but the plot in output is wrong. I guess I'm not accessing the solutions correctly... Jun 19, 2021 at 8:45
• That should be fine. Just make sure to plot with t too. Also, the default plot uses the interpolation to generate more values, so either use saveat directly or interpolate to get a denser plot. Jun 19, 2021 at 12:08
• How can I specify what t is? I see that the plot generated now goes from 0 to ~40 instead of 0 to 2. Did I lose time somehow when I used append? Jun 19, 2021 at 14:38
• How are you generating the data? soldadp(0.0:0.1:soldadp.t)` etc.? Jun 19, 2021 at 18:54
• You're on the right track. And something that may be helpful is the RecursiveArrayTools.jl DiffEqArray for having the t paired with the u for plotting. Jun 19, 2021 at 21:30