I think the insert function had a small bug too, e.g. the second cell would be sometimes set to zero instead of 50% of the parent.
This seems to work:
using OrdinaryDiffEq
function f(du,u,p,t)
for i in 1:length(u)
du[i] = p.α*u[i]
end
end
function condition(u,t,integrator)
return 1.0 - maximum(u) # has any cell hit 1?
end
function affect!(integrator)
u = integrator.u
idx = findall(y-> y >= 1.0, u) # find all cells which have hit concentration == 1
resize!(integrator,length(u)+length(idx)) # resize solution to accomodate all new divisions
for id in reverse(idx) # for each cell hitting 1, divide:
u[id+2:end] .= @view u[id+1:end-1]
u[id] *= 0.5
u[id+1] = u[id]
end
nothing
end
callback = ContinuousCallback(condition,affect!, rootfind=SciMLBase.RightRootFind)
u0 = [0.2]
tspan = (0.0,10.0)
p = (; α = 0.3)
prob = ODEProblem(f,u0,tspan, p)
sol = solve(prob,Tsit5(),callback=callback);