Optimising Julia ODE model

I'm currently trying to learn Julia but translating some of my code. Below, I exemplify with what kind of problems I'm usually confronted, i.e. small- to medium-sized systems of ODEs tested with a large battery of inputs and forcing functions. Here, I'm using ModelingToolkit.jl but I'm not sure whether this is the best way.

I'd appreciate any input on how to improve this code.

using DataFrames, DifferentialEquations, ModelingToolkit

# ModelingToolkit
@variables t xc(t) xp(t) u(t);
@parameters kel k12 k21;
D = Differential(t);

# Parameters
p = [kel => 1,
     k12 => 2,
     k21 => 1];

# Initial conditions
x0 = [xc => 0,
      xp => 0];

tspan = (0, 10);

# ODEs
odes = [D(xc) ~ -kel*xc - k12*xc + u + k21*xp,
        D(xp) ~           k12*xc + u - k21*xp];

# Stimulus function
# My function
function myfun(t, dose, duration)
    if t >= 0 && t <= duration
        dose
    else
        0
    end
end
@register myfun(t, dose, duration)

stim_size = [1, 2, 3];
stim_length = [1, 2, 3];

# Loop over stimulations
sol = DataFrame();

for i in 1:length(stim_size)
      for j in 1:length(stim_length)

            # Define ODE
            @named sys_ij = ODESystem([u ~ myfun(t, stim_size[i], stim_length[j]); odes])
            prob_ij = ODEProblem(structural_simplify(sys_ij), x0, tspan, p)
            sol_ij = solve(prob_ij)
      
            # Convert to dataframe and add stimulations
            sol_ij_df = DataFrame(sol_ij)
            sol_ij_df[!, :size] = repeat([stim_size[i]], size(sol_ij_df, 1))
            sol_ij_df[!, :length] = repeat([stim_length[j]], size(sol_ij_df, 1))

            # Combine
            append!(sol, sol_ij_df);
      end
end