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))
```