We present a rigorous approach, based on the concept of continuous thermomajorization, to algorithmically characterize the full set of energy occupations of a quantum system accessible from a given initial state through weak interactions with a heat bath. The algorithm can be deployed to solve complex optimization problems in out-of-equilibrium setups and it returns explicit elementary control sequences realizing optimal transformations. We illustrate this by finding optimal protocols in the context of cooling, work extraction, and catalysis. The same tools also allow one to quantitatively assess the role played by memory effects in the performance of thermodynamic protocols. We obtained exhaustive solutions on a laptop machine for systems with dimension d≤7, but with heuristic methods one could access much higher d.