Thank you for your response!
I’ve considered forcing the system to solve directly for B by instantiating it with the desired values
@variables B(t)[1:n] = rand(n)
and then setting the gradient to 0.
[D(B[i]) ~ 0 for i in eachindex(B)]...
This does solve the indexing time problem (at the cost of slightly higher solving time), but I was hoping for a cleaner solution.
Based on this thread it seems there is some way to mark variables you want to track, but the closest I’ve found is the irreducible metadata tag mentioned here:
@variables A(t)[1:n] = zeros(n) B(t)[1:n] [irreducible = true]
Sadly it slows down the solver immensely, making it an even worse solution for my problem.