Solving via Continuation: Perfomance in a loop

[ns-ku]

Hello!

I am having a problem with continuation with BoundaryValueDiffEq.
To solve two-point BVPs by continuation, I increase a parameter $C$ in a for-loop (code below).
A previously obtained solution is used as a new guess.
As was in the answer for a similar question here, I take advantage of DiffEqArray(sol.u, sol.t) to make it.

Here is the code I have (based on the pendulum code given as an example):

using BoundaryValueDiffEq, DifferentialEquations
using Plots

function main()

    tspan = (0.0, 1/2)
    C = 0.5
    g = 2
    p = (C, g)

    function test!(du, u, p, t)
        (C, g)  = p
        θ = u[1]
        dθ = u[2]
        du[1] = dθ
        du[2] = -(C*θ^2/g + dθ + g*θ + C*cos(C*θ))
    end

    function bc2a!(resid_a, u_a, p)
        resid_a[1] = u_a[1] 
    end

    function bc2b!(resid_b, u_b, p)
        resid_b[1] = u_b[1] - 1/2
    end

    function ig!(t)
        a = 1
        u1 = a*t^2
        u2 = 2*a*t
        return [u1, u2]
    end

    # first solution 
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), ig!, tspan, p,bcresid_prototype = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)

    # the parameter variation is set here   
    a = range(1, 4, 10)

    # loop for continuation
    for ii ∈ a
        p = (ii*C, g)
        bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
        sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
        println(ii)
    end

end

main()

However, even for a relatively simple system, the solver gets stuck/extremely slow at $C = 3$.
If do the same but manually, without main() function, everything works fast and well:

using BoundaryValueDiffEq, DifferentialEquations
using Plots
tspan = (0.0, 1/2)
C = 0.5
g = 2
p = (C, g)

function test!(du, u, p, t)
    (C, g)  = p
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(C*θ^2/g + dθ + g*θ + C*cos(C*θ))
end

function bc2a!(resid_a, u_a, p)
    resid_a[1] = u_a[1] 
end

function bc2b!(resid_b, u_b, p)
    resid_b[1] = u_b[1] - 1/2
end

function ig!(t)
    a = 1
    u1 = a*t^2
    u2 = 2*a*t
    return [u1, u2]
end

# first solution 
bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), ig!, tspan, p,bcresid_prototype = (zeros(1), zeros(1)))
sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)

# the parameter variation is set here   
a = range(1, 4, 10)

# loop for continuation
    p = (a[1]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot(sol)
    p = (a[2]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol)
    p = (a[3]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol)
    p = (a[4]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol)
    p = (a[5]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol)
    p = (a[6]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol)
    p = (a[7]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol)
    p = (a[8]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol)
    p = (a[9]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol)
    p = (a[10]*C, g)
    bvp = TwoPointBVProblem(test!, (bc2a!, bc2b!), DiffEqArray(sol.u, sol.t), tspan, p, bcresid_prototype  = (zeros(1), zeros(1)))
    sol = solve(bvp, MIRK4(), dt = 0.02, adaptive = true, abstol = 1e-9)
    plot!(sol) 

I am new to Julia. Am I missing something on the syntax part?
Thank you very much for any help with this!

--
I'm using
BoundaryValueDiffEq v5.8.0
DifferentialEquations v7.13.0

[ChrisRackauckas]

The loop creates a new scope: is it actually updating the initial condition?