using SymbolicIndexingInterface as SII
sol[SII.all_variable_symbols(sol)]
would generate the full time series of all unknowns + observed. If you solved with saveat it should be evenly spaced already. Just DataFrame and that and use SII.all_variable_symbols to set the header.
We can give that a nicer utility function sometime in the future but that should just be a few lines of code for what you need for now.