Solving chain of ODE in Julia

Question

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

function param_dadp!(s2,d,a0,ϕ₀)
    #ϕ₀ = abs(find_time!(s2,d))
    #ϕ₀=2.0
    function def_dadp!(dv,v,p,ϕ)
        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,ϕ₀)
    probdadp = ODEProblem(def_dadp!,α₀,tspan,(s2,d))
    soldadp = solve(probdadp,Tsit5(),callback=cb)
end
plot(param_dadp!(81.0,-0.0009,0.083163,2.0))

# 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

function param_dadp!(s2,d,a0,ϕ₀)
    function def_dadp!(dv,v,p,ϕ)
        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,ϕ₀)
    probdadp = ODEProblem(def_dadp!,α₀,tspan,(s2,d))
    soldadp = solve(probdadp,Tsit5(),callback=cb)

    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
    init = [last(soldadp);1.0]
    classspan = (last(soldadp.t),ϕ₀)
    probclass = ODEProblem(classic!,init,classspan)
    solclass = solve(probclass,Tsit5())

    solu = append!(soldadp[1,:],solclass[1,:])
    solt = append!(soldadp.t,solclass.t)
    sol = DiffEqArray(solu,solt)
end
plot(param_dadp!(81.0,-0.0009,0.083163,2.0))
```

Answer 1

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".