I have a stiff system of 21 ODEs, for which I want to use callback functions on some time-dependent parameters.

Looking for some help on this, as I’m running into some error messages and I’m new to callbacks. Not entirely sure how to incorporate the piecewise initial conditions.

Basically, I’m looking for:

```
v1 = 0 and v2 = 0 if t < a + d
v1 = gh and v2 = g(1-h) if t >= a+d
ydot = zeros(21,1) if t < a
ydot = ydot if t >= a
```

This is my code:

```
using LinearAlgebra, DifferentialEquations, MATLABDiffEq
a = 25;
d = 5;
T = 2323;
V = 1;
c = 0.1;
d = 0.35;
e = 0.5;
f = 0.2;
g = 1.0/3.0;
h = 0.3;
k = 0.12;
l = 1.0/11;
m = 0.6;
n = 0.1;
p = n*l;
q = (1-n)*l;
z = 0.01;
s = 0.2;
t = q;
u = 0.5*l;
v = 0.592;
w = v*u;
x = (1-v)*u;
y = 0.23245;
z = 13;
function RHS(ydot,y,par,t)
ydot = 0.0*zeros(21,1);
y = 0.0*zeros(21,1);
par = 0.0*zeros(23,1);
v1 = 0.0;
v2 = 0.0;
a = par[1];
d = par[2];
T = par[3];
d = par[4];
c = par[5];
e = par[6];
f = par[7];
g = par[8];
h = par[9];
k = par[10];
l = par[11];
m = par[12];
n = par[13];
p = par[14];
q = par[15];
z = par[16];
s = par[17];
t = par[18];
u = par[19];
w = par[20];
x = par[21];
y = par[22];
z = par[23];
fM = d * y[1] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;
h = c * d * y[2] * (y[19] + e * (y[5] + y[8] + y[11] + y[14]) + f * y[17])/T;
ydot[1] = -fM - (v1*y[1]-v2*y[2]);
ydot[2] = -h + (v1*y[1]-v2*y[2]);
ydot[3] = fM - k*y[3]-(v1*y[3]-v2*y[4]);
ydot[4] = h - k*y[4]+(v1*y[3]-v2*y[4]);
ydot[5] = k*(y[3]-y[5])-z*y[5]-(v1*y[5]-v2*y[6]);
ydot[6] = k*(y[4]-y[6])-z*y[6]+(v1*y[5]-v2*y[6]);
ydot[7] = z*(y[5]+y[6])-k*y[7];
ydot[8] = k*(y[5]-y[8])-z*y[8]-(v1*y[8]-v2*y[9]);
ydot[9] = k*(y[6]-y[9])-z*y[9]+(v1*y[8]-v2*y[9]);
ydot[10] = z*(y[8]+y[9])-k*y[10] + k*y[7];
ydot[11] = k*(y[8]-y[11])-z*y[11]-(v1*y[11]-v2*y[12]);
ydot[12] = k*(y[9]-y[12])-z*y[12]+(v1*y[11]-v2*y[12]);
ydot[13] = z*(y[11]+y[12])-k*y[13] + k*y[10];
ydot[14] = k*(y[11]-y[14])-z*y[14]-(v1*y[14]-v2*y[15]);
ydot[15] = k*(y[12]-y[15])-z*y[15]+(v1*y[14]-v2*y[15]);
ydot[16] = z*(y[14]+y[15])-k*y[16] + k*y[13];
ydot[17] = m*k*y[14]-z*y[17]-(v1*y[17]-v2*y[18])-t*y[17];
ydot[18] = m*k*y[15]-z*y[18]+(v1*y[17]-v2*y[18])-t*y[18];
ydot[19] = (1-m)*k*y[14]-(z+s)*y[19]-(v1*y[19]-v2*y[20])-(p+q)*y[19];
ydot[20] = (1-m)*k*y[15]-(z+s)*y[20]+(v1*y[19]-v2*y[20])-(p+q)*y[20];
ydot[21] = (1-m)*k*(y[14]+y[15]);
return ydot
end
function condition(out,y,t,integrator) # Event when event_f(u,t) == 0
out[1] = integrator.t - integrator.par[1] - integrator.par[2];
out[2] = integrator.par[1] - integrator.t;
end
function affect!(integrator, idx)
if idx == 1
v1 = integrator.par[8]*integrator.par[9];
v2 = integrator.par[8]*(1-integrator.par[9]);
else
ydot = 0.0*zeros(21,1);
end
end
cb = VectorContinuousCallback(condition,affect!,2);
y0 = [2323 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1];
tspan = (0.0, 100.0);
prob = ODEProblem(RHS,y0,tspan);
sol = solve(prob, MATLABDiffEq.ode15s(), callback = cb, saveat = 1.0)
```

and the error with stack trace:

```
ERROR: LoadError: BoundsError
Stacktrace:
[1] getindex at ./number.jl:78 [inlined]
[2] RHS(::Array{Any,2}, ::Array{ModelingToolkit.Operation,2}, ::Array{Any,1}, ::ModelingToolkit.Operation) at /Users/dc/desktop/is/post.jl:60
[3] (::ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{Any,2}, ::Vararg{Any,N} where N) at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/diffeqfunction.jl:248
[4] modelingtoolkitize(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}) at /Users/dc/.julia/packages/ModelingToolkit/wy693/src/systems/diffeqs/modelingtoolkitize.jl:19
[5] __solve(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::MATLABDiffEq.ode15s, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}; saveat::Float64, timeseries_errors::Bool, reltol::Float64, abstol::Float64, kwargs::Base.Iterators.Pairs{Symbol,VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64},Tuple{Symbol},NamedTuple{(:callback,),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}}}}) at /Users/dc/.julia/packages/MATLABDiffEq/b8nHE/src/MATLABDiffEq.jl:47
[6] solve_call(::ODEProblem{Array{Int64,2},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(RHS),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::MATLABDiffEq.ode15s; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:callback, :saveat),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64},Float64}}}) at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:96
[7] #solve_up#454 at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:122 [inlined]
[8] #solve#453 at /Users/dc/.julia/packages/DiffEqBase/1dZDc/src/solve.jl:106 [inlined]
[9] top-level scope at /Users/dc/desktop/is/post.jl:106
[10] include(::Module, ::String) at ./Base.jl:377
[11] exec_options(::Base.JLOptions) at ./client.jl:288
[12] _start() at ./client.jl:484
in expression starting at /Users/dc/desktop/is/post.jl:106
```

Thanks for any help!