我正在解决两个不同的 ODE，它们的解决方案需要匹配。我目前正在手动执行此操作，效果很好，但我想自动化此过程。第二个 ODE 从第一个解的最后一个值中获取一个初始条件和初始时间。我在下面放了一个MWE。现在，我天真地尝试将第一个和第二个 ODE 放在一起，如下面的第二个代码所示（给出 StackOverflowError）。不过，这似乎不是要走的路。有人知道该怎么做吗？

编辑：第二个代码现在是解决方案

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