We report predictions for the suppression and elliptic flow of the $Upsilon(1S)$, $Upsilon(2S)$, and $Upsilon(3S)$ as a function of centrality and transverse momentum in ultra-relativistic heavy-ion collisions. We obtain our predictions by numerically solving a Lindblad equation for the evolution of the heavy-quarkonium reduced density matrix derived using potential nonrelativistic QCD and the formalism of open quantum systems. To numerically solve the Lindblad equation, we make use of a stochastic unraveling called the quantum trajectories algorithm. This unraveling allows us to solve the Lindblad evolution equation efficiently on large lattices with no angular momentum cutoff. The resulting evolution describes the full 3D quantum and non-abelian evolution of the reduced density matrix for bottomonium states. We expand upon our previous work by treating differential observables and elliptic flow; this is made possible by a newly implemented Monte-Carlo sampling of physical trajectories. Our final results are compared to experimental data collected in $sqrt{s_{NN}} = 5.02$ TeV Pb-Pb collisions by the ALICE, ATLAS, and CMS collaborations.