Discrete sub-solutions for DEs

[YingboMa]

MTK models with synchronization and discrete subsystems need to store discrete solutions to be able to reconstruct the entire dynamics. It's easy to add a subsolutions field in ODESolution, but how should we attach the sub-solutions after __solve? Currently, MTK generates code that pushes the discrete states into a saving vector inside periodic callbacks. Should we add a higher order callback function like solution_post_processing after __solve in DiffEqBase to attach the sub-solutions?