In neutrino oscillations, a neutrino created with one flavor can be later detected with a different flavor, with some probability. In general, the probability is computed exactly by diagonalizing the Hamiltonian operator that describes the physical system and that drives the oscillations. Here we use an alternative method developed by Ohlsson & Snellman to compute exact oscillation probabilities, that bypasses diagonalization, and that produces expressions for the probabilities that are straightforward to implement. The method employs expansions of quantum operators in terms of SU(2) and SU(3) matrices. We implement the method in the code NuOscProbExact, which we make publicly available. It can be applied to any closed system of two or three neutrino flavors described by an arbitrary time-independent Hamiltonian. This includes, but is not limited to, oscillations in vacuum, in matter of constant density, with non-standard matter interactions, and in a Lorentz-violating background.