This is all a bit confusing, since the differential equations are not them selves functions of ‘C1’ and ‘C2’, so a differential equation solver might not be necessary. Jus multiply ther RHS of those equations by the current value of ‘t’ an be done with it.
With respect to ‘r’ since it is a function of ‘C1’ and ‘C2’ that have initial conditions, calculate it, ‘P1’ and ‘P2’ first in each loop iteration.
At least, knowing only what is in the problem statement, that’s what I would do.
.
