I am trying to solve 4 coupled differential equations the first two of them should always be positive while the 3rd and 4th can be positive or negative. I found that I can impose a positive domain on all four variables using
sol = solve(prob,AutoVern7(Rodas5()),callback=PositiveDomain(),reltol=1e-8,abstol=1e-8)
But is there a way to restrict the positive domain to a list of variables? Ideally, my plan is to work with a huge system where every third and fourth variable out of a few hundred has to be positive while others can be negative.